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CHAPTER I 


INTRODUCTION 

The effects of atmospheric turbulence on many of 
the modern sophisticated technological systems have become 
an important design parameter from both structural and per- 
formance aspects. Techniques for simulating atmospheric 
turbulence have therefore been developed in an attempt to 
provide reliable design criteria. Turbulence simulation is 
achieved by the generation of an analog or digital signal 
which has equivalent statistical characteristics to the 
true atmospheric turbulence. The degree of complexity and 
mathematical involvement of the models has increased with 
each successive generation, from a simple, one-component 
wind speed having a Gaussian distribution and a Dryden 
turbulence spectra to multi-component models having non- 
Gaussian probability distribution, more complex spectra, 
and simulation of other statistical properties, such as 
coherence. The simulated turbulence can then be used to 
predict the behavior of airplanes, bridges, buildings, etc., 
under the influence of a turbulent atmosphere. The purpose 
of this study is to investigate the effects of three turbu- 
lence models of various complexity on the performance of a 
simulated aircraft landing under an automatic controlled 
system. The work compares the aircraft landing under a 
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simple "Z-transform" simulation technique [1] , 1 a non- 
Gaussian simulation technique [2] , and a simulation which 
incorporates vertical coherence 13] . 

The basic simulation procedure is shown in 
Figure 1-1. Random signals are computer generated and 
passed through shaping filters to provide the output which 
is the simulated time history of turbulence. The degree of 
complexity of the system is internal to the filter system. 
Through the appropriate design of the filter, however, the 
turbulence signal output is designed to have certain 
statistical properties which are the same as those that 
have been measured and verified for atmospheric turbulence. 



Output 


Figure 1-1. Turbulence simulation system. 


This report first reviews the statistical 
properties of atmospheric turbulence, in particular the 
probability distribution, the spectra, and the coherence 
in Chapter II. The three different simulation techniques 
investigated in the study are then described in Chapter III, 

^Numbers in brackets refer to similarly numbered 
references in the Bibliography. 





and appropriate statistical analyses are carried out to 
verify validity of the models in Chapter IV. Finally, in 
Chapter V, the models are incorporated into a computer 
model of aircraft flight dynamics; and statistical landing 
results for simulated flights of an aircraft having charac 
teristics of a DC-8 are made for the different turbulence 
simulation techniques. The significance of the various 
degrees of sophistication introduced into the simulation 
technique on the landing performance of the aircraft is 


discussed . 


CHAPTER II 


ATMOSPHERIC TURBULENCE AND MODELS 

This chapter presents the basic expressions 
involved in mathematically describing atmospheric turbu- 
lence. In particular, the statistical properties such as 
the probability density distribution of wind speeds, the 
turbulence energy spectrum functions ,. and the spatial 
coherence are discussed. 

A. Atmospheric Turbulence Spectra 

In developing the atmospheric turbulence statistics, 
the wind velocity is written as 

U. = U. + u. , (2-1) 

111' 

where U^ is the instantaneous velocity, U. is the average 
velocity, and u represents the fluctuating component. The 
subscript i denotes the i direction of the wind velocity. 
The intensity of the fluctuations in wind speed is 
expressed by the root-mean-square value 


a . 

i 


= (up 


1/2 


( 2 - 2 ) 


The mean turbulence kinetic energy, TKE , is pro- 

2 

portional to the value of u' and includes contributions 
from all frequencies of eddies making up the turbulent 
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motion. If the time history of the velocity fluctuations 
is processed through a filter which passes only a small 
selected band of frequencies, Aoo, the mean square value of 
the wind fluctuation in that frequency range will be pro- 
portional to the turbulence kinetic energy of gusts or 
perturbation in the wind having that frequency. The value 
of the turbulence kinetic energy density function at the 
midpoint of the band is then given by 

E. .(to) = \T. uT / A to . (2-3) 

13 1 3 


The function E^ . (to) forms the spectrum tensor 
having nine components. The function E^(to) defined by 
Equation (2-3) for i = j is called the one-dimensional 
energy spectrum function. It follows that 


uT = / E . . (to) dco . 
x 11 


(2-4) 


Frequently the energy spectral density is normalized as 


^(u) = E ±i (to) / U 2 , 


(2-5) 


and thus , 


/ <j>. ( 10 ) do) = u? / U 2 = a? / U 2 
o 


2 

where a. is called the variance. 
1 


( 2 - 6 ) 
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For i ? j , the function (w) is related to the 
cross-correlations and is called the cross -spectrum. 

The energy spectrum tensor can be considerably 
simplified if the atmospheric turbulence is homogeneous and 
isotropic, for which there is no mean rate of transfer of 
momentum across shearing surfaces or, more specifically, 
the Reynolds or eddy shearing stresses, -puTuT , vanish. 

For this case, a three-dimensional energy spectrum function 
E (K) is defined as 

E (K) = 21TK 2 IE U (K) + E 22 (K) + E 33 (K) ] , (2-7) 

where K is the wave number defined as K = w/V, and V is the 
reference velocity? therefore, E^ (K) - VE^ (to) . 

Figure 2-1 shows the typical form of the three- 
dimensional energy spectrum. 



wavenumber 

K 


Figure 2-1 . Typical three-dimensional energy 
spectrum function [4] . 
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Physical argument coupled with dimensional analysis 
has provided basic insight into the behavior of E (K) . The 
mathematical form of E(K) adopted by von Karman [5] is 


E(K) 


( K / K q ) 4 

2 17/6 ' 
[ 1 + ( K / K q ) J 


(2-8) 


where c = (55/9) (A/p) a 2 , K = 1/(1. 339A), and A is the 
integral length scale defined by the correlation function R 
as 


A = / R(t)dt . (2-9) 

o 


4 

This von Karman spectrum function exhibits a K and 
K~5/3 j-, e havior a t i ow an £ high frequencies, respectively, 
which represents the limiting conditions of the spectrum 
for these wave number ranges. Also, the one-dimensional 
spectrum function, $ , becomes 


^(K) 


4> 2 (k) 


<j> 3 (k) 


2 

= a l 

2A i 

1 

7T 

2 5/ ^ 6 

[ 1 + ( 1.339 A ± K ) ] 

- a 2 

A 2 

1 + 8/3 ( 1.339 A 2 K ) 2 

" a 2 

TT 

T~IT76 

[ 1 + ( 1.339 A 2 K ) ] 

2 

— fX 

^3 

1 + 8/3 ( 1.339 A 3 K ) 2 

a 3 

TT 

, U“/6 

[ 1 + ( 1.339 A 3 K ) ] 


( 2 - 10 ) 
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where the subscripts 1, 2, and 3 represent longitudinal, 
lateral, and vertical fluctuations, respectively. 

Also, Dryden [5] points out that for laboratory- 
type flow the shape of the velocity correlation curve can 
be approximated by an exponential function. In this case 
the one-dimensional Dryden spectra become 


4> X (K) 


2 (K) 


<f> 3 (K) 


2A X 

1 

t r 

2 2 

i + 

^2 _ 

2 2 

1+3 A^T 

7T 

( 1 + A|k2 ) 2 

A 3 _ 

1+3 A 2 K 2 

TT 

( 1 + a 2 k 2 ) 


( 2 - 11 ) 


Simulation of atmospheric, turbulence usually 

-2 

employs the Dryden spectrum which exhibits a K fall-off 
at high frequencies and is, therefore, easier to handle 
mathematically, as described later. However, as Figure 2-2 
shows, the von Karman spectrum gives a better description 
of the measured data than the Dryden spectrum. 

In a recent paper describing the spectral 
properties of atmospheric turbulence over a flat homo- 
geneous field site, Kaimal [6] shows that with proper 
nondimens iona li zation , spectra in the stable atmospheric 




-4 -3-2-10 1 2 

10 10 10 10 10 10 10 1 ( 


0.55 K/k 
m 


Figure 2-2. Longitudinal velocity spectra compared with 
the von Karman and Dryden spectra (wave- 
number normalized by the peak scale K of 
the vertical spectra) [8], m 
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surface layer can be reduced to a universal curve having 
the empirical formula 


. 0.164 (n/n) 

— <P (f ) = K'/'k ' (2-12) 

a 1 + 0.164 ( n/o 0 ) / 

where f is the cyclic frequency, n = fz/V is the reduced 
frequency, z is height, and r, Q is a scaling parameter 
related to the atmospheric stability condition. Figure 2-3 
shows n Q plotted both against the ratio of height to Monin 
Obukhov length scale, z/L, and against the gradient 
Richardson number, Ri. For neutral conditions, the values 
of n Q recommended by Frost, et al. [7] are 

n n = 0.0144 
U 1 

ri n = 0.0265 (2-13) 

U 2 

=0.0962 . 

U 3 

Figure 2-4 shows experimental data for atmospheric 
turbulence compared with Equation (2-12) . 

B. Filter Theory for Simulation Applications 

The principal use of power spectral density 
functions is to establish the frequency composition of the 
data which, in turn, bears an important relationship to the 
basic characteristics of a physical system exposed to inter- 
action with the turbulence. For example, consider a system 





13 


exposed to interaction with the turbulence. For example, 
consider a system as shown in Figure 1-1, page 2, with a 
frequency response function H(a>). Assume that a stationary 
random signal with a power spectral density function <J> x ((o) 
is applied as an input to this system. The output from the 
system will be a stationary random signal with a power 
spectral density function <£y((o) given by 


<P y M 


H (to) 


2 <!> x (u) • 


(2-14) 


One can obviously design a filter function to 
obtain the form of the output spectrum <j> (w) desired for 
any given known input (p (to) . For example, if the desired 
output is the longitudinal or lateral component of the 
Dryden spectrum given by Equation (2-11) , then Equation 
(2-14) becomes 


4> x (w) H 1 ((0) 


° 1 2A 1 1 

VlT 1 + ( AjU/V ) 2 


, a 2 A, 1 + 3 ( A-id/V ) 2 

♦x ( “) -~Vir H 

1 1 + ( a 2 w/v r ] 


(2-15) 


Introducing the input as wide-band white noise for 
which the spectrum 4> x (co) is a constant, adjustable to unity, 
gives 


14 


where 


H 1 (a)) | 2 = 


H 2 (w) 


2 2 
+ a) 


2 2 

c 2 ( b + go ) 
2 

, 2 ^ 2 “ 

( a 2 + a) ) 


c-, = 


o 2 

2a i p i 

7T 


a-, = 


V 

A, 


(2-16) 


o 2 
3 a 2°2 
7 T 


b = 


/3 


a,-, = 


V 


2 A. 


One solution for H(o>) is 


✓cj 

H (03) = -■ ■ 

1 a., + 303 


H 2 (oj) = 


/cj ( b + j o) ) 
[ a 2 + j 03 ] 2 


(2-17) 


The complex polar notation gives the frequency response 
function in terms of a gain factor | H ( 03 ) | and a phase 
factor 0 ( 03 ) as 


= |H( U ) |e- je( “ ) 


H(to) 


(2-18) 
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Thus t 


H-^w) = 




r 2 , 2 , 

[ a^ + to ] 


I7* 


e x (to) = tan ^ ( w/a 1 ) 


and 


H 2 (w) = 


, ,2 , 2 » 

( b +oo ) 

, 2 2 . 2 
( a 2 + w ) 


^/2 


(2-19) 


0 2 ( oo ) = tan 


-1 


( 2a 2 b - ( a 2 - m 2 ) ) 
2a 2 w 2 - ( a 2 - to 2 ) b 


These gain factors and phase factors are plotted in 
Figure 2-5. Because of the irrational form of the Kaimal 
and von Karrnan spectra, approximate forms of H(tu) are used. 
Equation (2-20) is a modified Kaimal spectrum with a Dryden 
form: 


f <!>(£) _ 
2 


= 1.9 


( n/n 0 ) 


(3.8) + 


( n/n Q ) 


( 2 - 20 ) 


Figure 2-6 shows a comparison of this approximate form with 
the original spectrum. 







f<j> (f )/o 2 



100.0 


The Kaimal spectrum. Equation (2-12) , 
compared with the approximation form, 
Equation (2-20) . 


From Kaimal [6] a definition of A is given as 
0.041 (z/ri Q ) , which gives 


n/n o 0.041 V • 


Substituting Equation (2-21) into Equation (2-20) and 
replacing f with o)/2tt gives 


♦<“) = -2 ~ 2 
a +o) 


where 


0.156 a irV 


0.312 • ttV 
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Recalling Equation (2-16) , it follows that the frequency 
response function H(w) for the approximate Kaimal spectrum 
can be solved in the same manner as that used to obtain 
Equation (2-17) . 

The approximate von Karman spectrum, as shown in 

19] is 

1.942 a 2 A I 1 + ( 3.496 A f ) 2 ] 

<H f) = ^ \ , (2-23) 

[ 1 + ( 3.496 A f ) 2 ] 

which can also be written as 

<f>(w) = c ( b - +ai — L . (2-24) 

[a 2 + w 2 ] 


The filter function H(to) in this case is therefore similar 
to the lateral component of Equation (2-17) . Figure 2-7 
compares the approximate expression with the actual 
von Karman longitudinal spectrum given by Equation (2-10) . 

C, Multi-Filter System for Non-Gaussian Turbulence 
and Turbulence with Interlevel Coherence 

Additional turbulence characteristics can be taken 
into account with more complex filter systems. A filter 
system which gives a non-Gaussian turbulence output and a 
turbulence output with interlevel coherence, respectively. 


is described in this section. Both simulations use the 
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same Dryden spectrum. The frequency response function for 
each filter can no longer be derived from the simple 
relation given in Equation (2-14) . 

The non-Gaussian output filter. Figure 2-8, com- 
bines three filter functions, H , H, , and H , for each of 
the three wind speed components being simulated. This 
nonlinear system allows the probability density function to 
be adjusted through the parameter r. Figures 2-9 and 2-10 
show the comparison of measured atmospheric turbulence with 
a Gaussian distribution curve, and the comparison of the 
same turbulence with the probability density function simu- 
lated by selecting different values of r, respectively. 


In ? ut Filter °? tpu * 

noise signal 


Gain 


Simulated 

turbulence 



g (t) 


Figure 2-8. The non-Gaussian multi-filter system. 


To determine the frequency response functions, 

H , H, , and H , the correlation function is first derived 
a 9 b ' c 

from the Fourier transform of the spectrum to be simulated: 
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R(t) = / 4> (K) exp (iVKt) dK . (2-25) 

— CO 

Considering the Dryden spectra. Equation (2-11) , 
the correlation becomes 


R-^t) = exp [ - V/A 1 1 1 J 


R 2 (t) = o 2 2 [ 1 - ] exp [ - V/A 1 1 1 ] . 


(2-26) 


The correlation function of the simulated turbu- 
lence time history can also be determined from (see 
Figure 2 r- 8 ) 


R (t) = Elg(t) g (t + t) ] , (2-27) 

y y 


where 


g (t) = a (t) b (t) 


[ 1 +r ] 


1/2 


+ c(t) 


[ 1 + r J 


1/2 


(2-28) 


Since a(t), b(t), and c(t) are independent, random 
processes with zero mean, the correlation function can be 
written as 


R gg (t) =R aa (t) *bb (t) 


[ 1 + r ] 


o-+ R 

2 , cc 


T~ ' 

[ 1 + r ] 


(2-29) 


and each function can be related to a spectrum function 
through Equation (2-30) : 
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00 

R aa (t) = / 4> aa (K) exp (iVKt) dK , (2-30) 


where $ (K) can be replaced by the filter function H(iVK), 

da 

i .e . , 


|H a (iVK) | 2 = <t> aa (K) 


(2-31) 


Similar results apply to and R cc - Reeves, et al. [2] 

assume the general form of the response function to be 


N 1 

H a (s) = I~TDTi- 


H b (s) = 


n 2 + n 3 s 
( 1 + d 2 s t 


H (s) = 
c 


N 4 + n 5 s 
( i + d 3 s Y 


(2-32) 


where s is the Laplace transform variable s = joj. 

Substituting Equation (2-32) into Equation (2-31) 
and transforming as shown by Equation (2-30) yields the 
following correlation functions: 
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R aa (t) " 2E[ exp [_ 


N, N_ - 2 

Kbbtt) = (jJ-) { |t| I < fjJ) - (i) ] 


+ [ &7 + D 2 ( ST 1 3 } exp ] (2-33) 


R cc (t) = { N 1 <i£> - <67 


N 


+ [ / + D 3 ( n 7 ) ]}ex P < 

3 5 


I tl 


These particular expressions for the response functions 
have been chosen because they will give the correct Dryden 
spectrum. Also, to eliminate the r appearing in Equation 
(2-29) , the following choices are made for the constants in 
Equation (2-32) : 


N, = 


Nr = 


4a l V ' N 2 = 1 *° ' N 3 = ' N 4 = a, (^) 


2A 
1 1 V 


1/2 


a l ( 


2A' 


1/2 


V’ 


D= 2A d= 2A d= A 

U 1 V ' 2 V ' 3 V 


The resulting correlation functions of the model become 


R aa (t) R bb (t > = R cc (t) 


R (t) 

gg 


2 r v 4 _ i 

exp [- J t ] 


(2-34) 
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This is the desired form of the longitudinal gust corre- 
lation function for a Dryden spectrum. Similarly, the 
lateral and vertical components can be simulated correctly 
by proper choice of constants. These constants are listed 
in Table 2-1, 

Atmospheric turbulence near the ground has been 
shown to have strong coherence between vertically separated 
layers of air. Many researchers [11, 12, 13] have found 
that the coherence near the ground (z < 150 m) behaves a§ 

Y o (K, Az) = exp [-a |KAz| ] . (2-35) 

Therefore, this statistical property should be included in 
a valid turbulence simulation. 

Perlmutter, Frost, and Fichtl [3] developed a 
multi-filter system to provide turbulence simulation 
including coherence function. They defined the two-point 
spectral function cj) (K , , z £ ) at different height z^,z^ 

to be 

o P * 

t(K,z,,z 2 ) =0(K,z,)H*(K,z,) l D m (K,z 1 )D in (K,z 2 ) , (2-36) 

m=-P 

where H(K,z) is the frequency response function, and D(K,z) 
is a level frequency factor, as shown in Figure 2-11. 
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Table 2-1. Transfer Functions for the 
Non-Gaussian Model 


H. 


Longitudinal 


1+2 (A) s 


1+2 4 ) s 


<7 (— ) 

1 ^ V J 


1/2 


l+C^s 


Lateral 


a 2 (128) 1/2 (A) 

1+ 2 (A)s 


(1+2^3) 2 


a 1 / 2 A 

° 2 <$ Cl+^js) 

(l+(A)s) 2 


a 3 (128) 1/2 (A) 

1+2 (A) s 


a+iAs ) 2 


a 3 (A) fl+/fe) 

(i+(A)s ) 2 


1/2 


Vertical 
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Figure 2-11. The interlevel coherence model. 


The coherence function y(K,z^,z 2 ) is defined as 


4) (K,z lf z 2 ) 4 )*(K,z 1 ,z 2 ) 
Y(K ' z 1 ' z 2 ) = 4 >(K,z 1 ,z 1 )((.(K,z 2 ,z 2 ) _ • 


(2-37) 


The level frequency factor is therefore assumed to have the 
form 


D (K , z) = c exp [ j k z d ] , 
m m c J m 


(2-38) 


where 


d = 


mu 


m (KAz). 


max 


Therefore, Equation (2-36) becomes 


<MK,z ,z ) =A 2 H(K,z )H*(K,z 2 ) l c m exp[jk(z 1 -z 2 )d m ] 
x z m=-P 
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Substituting this into Equation (2-37) yields 



Notice that c and d are equal for D and D y therefore, 
mm m — m 

Equation (2-39) can be rewritten as 


y (K, Az) 



(2-40) 


The error between this equation and the true wind 
coherence exponential function given by Equation (2-35) is 


E 


rr 



Y (K, Az) 


y Q (K,Az ) ] d£ , 


where 


£ 


KAz 


(KAz). 


max 


(2-41) 


Values of c which cause E to be a minimum are found by 
m rr 

setting dE rr /dc m = 


0 . 


The result is: 
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To determine the value of .the parameter A in 
Equation (2-36) , consider that Equation (2-36) must reduce 
to the. one -level auto-spectral function when Az = 0 . 
Equation (2-36) then takes on the new form of 

P 

<HK,z 1 ,z 1 ) = A^H(K,z^)H*(K,z^) [c^ + 2 l c^] . (2-43) 

m=l 

Comparing Equation (2-43) with a single filter system for 
which the spectrum is 

(HKjZ^ = |h(K,z 1 ) | 2 , (2-44) 


then A can be obtained as 


A 


2 



(2-45) 


Therefore, the two-point interlevel coherence model reduces, 
as necessary, to the one-level auto-spectrum at Az = 0. 



CHAPTER III 


DIGITAL SIMULATION TECHNIQUES 

This chapter describes the basic methods and 
problems associated with digital simulation of random 
processes. Topics discussed are the generation of a 
random signal and the calculation of the necessary trans- 
formation functions used in later work. 

A, Generation of a Random Signal 

The method used to generate random numbers with a 
digital computer is based on the formula 

x n+ ^ = a x n + b (mod m) . (3-1) 

This technique is called the mixed congruential 
method and provides a simple and fast computation pro- 
cedure. It also has the advantage that it can be repeated 
for several cycles by using different values of a and b. 
Moreover, it can be analyzed theoretically. The period at 
which the random numbers repeat cannot be greater than the 
module number m; thus, m is usually selected to be the 
largest integer possible with the capability of the com- 
puter used. The value of m in this study is 2,147,483,647, 

31 

which is equal to (2 - 1) , the largest possible number 

allowed on the IBM 360/65 system. 


31 



32 


Frequently, it is desired that the sequences 
generated be bounded by zero and unity and be uniformly 
distributed within this interval. This bounding require- 
ment is accomplished easily by normalizing the output data 
with respect to the largest possible number m. 

x* = x /m , (3-2) 

n n 7 ' v ' 


where is generated from Equation (3-1) . 
Since 


therefore , 


0 < x < m ; 

- n - 


i 

0 < x < 1 . 
n - 


The constants a and b in Equation (3-1) are selected to 
provide speed of computation and good statistical 
properties, such as uniform distribution and maximum period. 

Results presented by Chambers [14] and Hamming [15] 
show that a maximum period can be achieved if a and b 
are selected as 


fa =4*1+1 I = 1,2,3,... 
b = odd 


or 


(3-3) 


a = 8 * I ± 3 
b = 0 . 
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The computer program RANDU used for generating 
random numbers was combined in subroutine GAUSS, as dis- 
cussed in the following text. The value a is selected 
equal to 65539 which follows the second rule of Equation 
(3-3) by setting I equal to 8192 and provides the best 
statistical results using the IBM 360 computer. 

As mentioned earlier, Gaussian distributed numbers 
are required as the input white noise. The probability 
density function of the numbers must, therefore, have the 
form 


__ 2 

P(x) = — - — exp I- ] , (3-4) 

a /2¥ 2a 

x x 

where a is the standard deviation and x is the mean value 
x 

of the variate. 

If x|(i = 1,2,...N) are N uniformly distributed 
numbers over the interval 0 to 1 , the central limit theorem 
yields the following formula: 

12 V 2 N 

x\' = 1^3 [I X*. -N/2]+x , (3-5) 

3 x in i=1 

where x ( j = 1,2,...M) is a set of random numbers having an 
approximate Gaussian distribution with the mean value equal 
to x, and the standard deviation equal to a^. 

The value of N indicates the number of terms used 
for each output. Figure 3-1 compares the approximation 



Figure 3^1. The Gaussian distribution 

compared with the approximation 
Equation (3-5) [16] . 
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given by Equation (3-5) for 5 terms and 12 terms with the 
true Gaussian distribution function. 

The computer program for generating Gaussian random 
numbers also appears in the IBM scientific subroutine 
package. The N value of Equation (3-5) is selected equal 
to 12 in order to simplify the computations. Equation (3-5) 
then reduces to 

N _ 

Xi - a r I x; - 6 ] + x . (3-6) 

3 x i=l 1 

The computer subroutine is called GAUSS and is given in 
Appendix A, Section A-l. 

B, Calculation of the Required Discrete 
Transformations 

The filter system shown in Figure 3-2 consists of 
an input and an output history function x(t) , y(t) , 
respectively, and a system response function h(t). 


x(t) 


h (t) 


y (t) 


Figure 3-2. Simple filter system. 


The output y(t) is given by the convolution integral 

OO 

y (t) = f h(x) x(t-T) dx . 


— OO 


(3-7) 
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The Fourier transformation of Equation (3-7) is [17J 

Y(u) = H(u>) X(to) . (3-8) 

Since the filter function H (o> ) can be designed to simulate 
most forms of the spectrum desired, as described in 
Chapter II, Section B, the calculation of a simulated time 
history of turbulence relies mainly on the Fourier trans- 
formation of the input function x(t) to X(co) . The inverse 
transformation of the Y(w) back to y(t) is also of equal 
importance to the computational effort. These two pro- 
cedures will be discussed in detail in the following, as 
they pertain to transformation of digital data. 

The Discrete Fourier Transformation 

For cyclic frequency f, the infinite range Fourier 
transform of a real valued or a complex valued record x(t) 
is defined by 

CO B 

X (f ) = / x ( t) e~ 32irft dt . (3-9) 

*-oo 

It is physically impossible to calculate this transfor- 
mation since an infinite number of digitized data does not 
exist in reality. However, by restricting the limits to a 
finite time interval of x(t), say in the range (0,T), then 
the finite range Fourier transform will exist as defined by 
Bendat and Piersol [17] 
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X(f,T) - / x(t) e ^ 2irft dt . (3-10) 

0 

Assume now that. x(t) is sampled at N equally spaced 

points separated by a distance At apart f where At has been 

selected to produce a high cut-off frequency to avoid 

aliasing. The cut-off frequency (Nyquist frequency) and 

the aliasing problem are discussed in a later section. The 

discrete values of x thus become 

n 

x n = x (n At) n = 1,1,2, ...N-l , (3-11) 

and the discrete version of Equation (3-10) becomes 

X(f,t) = At j 0 X n ex P 2 * f nAt] . (3-12) 


A basic frequency f is defined as 

f = ~ , (3-13) 

o T f 

and the discrete frequency value f is defined in terms of 

f as 
o 

f = m f = m = 0,1,2, ...N-l . (3-14) 

m o NAt ' ' ' 


The discrete frequency f is substituted into the 
discrete Fourier components, X^, defined as 


X 


m 


X(f T) 
' m 


At 


(3-15) 
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Therefore , 


N-l 

X m = I x n ex P 3 
n=0 


. 2 7rmn 
N 


(3-16) 


To simplify the notation, let 


W = exo r- i 2 7111111 l 
mn exp 1 3 — J 


(3-17) 


Then Equation (3-16) becomes in tensor form. 


X 


m 


W 

mn 


x 


n 


(3-18) 


This is the equation for the discrete, finite- 
grange Fourier transform. No difficulty is involved in 
carrying out calculation of the Fourier transform with the 
equation using a digital computer. Equation (3-8) can now 
be written in discrete form as 


Y 


m 


H X 
m m 


(3-19) 


This equation gives the discrete Fourier components of the 
output to be simulated. To invert the output Fourier 
components into the desired time domain, the discrete 
inverse Fourier transform is applied, which can be derived 
similar to Equation (3-18) , The result is 
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where 


y = i W* Y 
J n N ran m 


W mn = -p [ j 2™» ] 


(3-20) 


These three steps complete the calculation; that 

is, the simulated discrete data y n are obtained by one 

Fourier transform, Equation (3-18) , one multiplication, 

Equation (3-19) , and one inverse Fourier transform. 

Equation (3-20) . To calculate Equation (3-18) or 

. 2 

Equation (3-20) requires N times of complex multiplication 

and N(N-l) complex additions to transform N discrete values. 

Since the value of W repeats for certain combi- 

mn 

nations of m and n, some additions and multiplications can 
be eliminated. The Fast Fourier Transform (FFT) is an 
algorithm based on this principle which allows computation 
of the transformation much more rapidly than the direct 
method in Equation (3-18) or Equation (3-20) . The FFT per- 
forms a series of computations by rearranging the matrix 

W to save the calculation of dual node pairs. For N 
mn 

discrete points, where N is a power of 2, namely, 

N = 2 P , (3-21) 

the FFT method requires only P»N/2 complex multiplications 
and P*N complex additions. Thus, the approximate ratio of 
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computing time for the direct method compared to the FFT 
method is 


Direcrtf Method _ 2 n 
FFT ~ P 


(3-22) 


Figure 3-3 illustrates the difference between the 
number of multiplications required for each method. 



Figure 3-3. Comparison of the number of multiplications 

required for a Fourier transform by the direct 
algorithm and the FFT algorithm. 


The Z-Transf ormation 

Another useful discrete transformation technique is 
the Z-transformation which is defined as 


= J x (m At) z 
m=0 


m 


f 


X(z) 


(3-23) 
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where X(z) is called the Z-transform function of X(t). 
Actually , the Z -transformation is the Laplace transform of 
digitized data. It can be derived from the Laplace trans- 
formation defined as 


X(s) = / x(t) e dt . 
o 

Consider the data x(t) digitized as 

oo 

x’ (t) = £ x (m At) 6 (t-m At) , 

m=0 

where 

6 (a) = 1 for a = 0 

= 0 for a ^ 0 . 

Taking the Laplace transform of x f (t) gives 

OO 

X' (s) = / x' (t) e~ st dt . 
o 


(3-24) 


(3-25) 


(3-26) 


Substituting Equation (3-25) into Equation (3-26) results 
in 

X'(s) = f x (m At) e -s(m At) . (3-27) 

m=0 

The Z-transform X(z) is equal to X'(s) if z is set equal 

Q A ^ 

to e , and Equation (3-27) becomes 
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X(z) = x'(s) = J x e~ m(s At) f (3-28) 

m-0 m 

where 

x m = x (m At) . 

The transformation is now defined as a function Z, 
X(z) - Z I x (t) ] . (3-29) 


Therefore, the transf ormation of x(t + At) becomes 


Z I x m+1 ] = I * m+ l z 


-m 


m=0 


= z ( J n X m+l Z 

m=0 


- (m+1) 


/V -m \ 

= z ( > x z - x ) 

m o 

m= 0 

= z ( Z [ x ] - x ) 
m o 


(3-30) 


in the same manner it can be shown that the Z -transfor- 
mation of { x , - x ] is equal to 
m+1 m 

Z [ x , - x ] = (z-1) Z [ x ] - z x . (3-31) 

m+1 m mo 

This equation can be used to determine x m+ ^ from x^ 
by the appropriate inverse transformation. The application 
of Equation (3-31) to simulate turbulence is as follows. 



43 


The basis of the Z-transform technique is to 

establish the difference equation for the output 

Therefore, if the Z-transform of the filter function h(t) 

in the system shown in Figure 3-2, page 35, can be found, 

then the output Y n+ -^ can be calculated directly from the y n 

and x values, 
n 

An example determination of the filter function 
h(t) for simulation of the Dryden spectrum is given in the 
following. 

The Fourier transform is (see Chapter II, Equation 

(2-17) ) 


H (a)) 


/c“ 

VC 1 


+ a. 


; 


which also can be expressed in Laplace form by setting 
s = j cu (see [18] ) 


H (s) 


/c~ 

VC 1 


s + a. 


(3-32) 


Now, consider a unit step function x(t) as the input to the 
system where 


x(t) =1 t > 0 

= 0 t < 0 . 


(3-33) 


The Laplace transform x(t) can be obtained from Equation 
(3-24) as 


X(s) = L [ x(t) ] = 1/s . 


(3-34) 
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For a linear system, the Laplace transform of the output is 
related to the input by 

Y (s) = H (s) * X(s) . (3-35) 


Therefore, the output Y(s) of passing a unit step input 
through a Dry den spectrum filter would be 


Y(s) 


s+a n 


1 

s 


(3-36) 


Partial fractions give 


Y(s) 



1 

s+a l 


) . 


(3-37) 


Taking the inverse Laplace transformation provides the 
output y(t) 


y(t) 


/c~ 

-1 V ^1 
L 1 [ Y(s) J - — ± 

a l 


(1 



(3-38) 


The Z-transform of the input x(t) and of the output 
y (t) can be obtained by substituting Equations (3-33) and 
(3-38) into Equation (3-23); hence. 


X(z) 


2 

z-1 


Y(z) = ( 

a l 


z 

z-1 


-a , At 
z-e J- 


(3-39) 
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Since 


H ( z) = Y ( z ) / X ( z) 


(3-40) 


the Z-transform of the Dryden spectrum filter function 
becomes 


H ( z) 



-a. At 
1-e 1 

-a 1 At 

z-e 


) . 


(3-41) 


For a general input X(z) the output of the filter system 
can be written as 


Y(z) 



-a., At 
1-e 1 

“ a l At 

z-e 


) X (z) . 


(3-42) 


Cross multiplication results in 

(z-A)Y(z) = BX(z) , 


(3-43) 


where 


-a, At /cT -a, At 

A = e , B = — ( 1-e ) . 

a i 


Now, applying Equation (3-31) and setting y(0) - 0, the 
left-hand part of Equation (3-43) becomes 


Z 1 [ y m+l " Ay m 3 = (Z ‘ A)Z [y m 3 

(3-44) 


(z-A)Y(z) . 
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Hence, Equation (3-43) becomes 

Z =BZ Ix m ] • < 3 - 45 > 

Taking the inverse Z-transf ormation , Equation (3-45) 
becomes 


and 


Vn " Ay m = Bx 


m 


m+1 


= Ay + Bx 


m 


m 


(3-46) 


This equation for calculating the discrete output 
y(t) is simple and very convenient for digitized data. 

Also, because the coefficients A and B are predetermined, 
this approach does not require calculation of the filter 
function, H, over all values of frequency in each compu- 
tational step, as required in the FFT technique. For these 
reasons, the Z-transformation method is even more efficient 
than the FFT method. Only 2N real multiplications and N 
real additions are required to simulate N discrete 
components of a given time history. If N equals 1024, for 
example, the FFT method requires ten times more multipli- 
cations than does the Z-transform technique. 
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C. General Consideration in Digitization and 
Discrete Transformation Calculation 

There are several required operations for the 
processing of random data. They are heavily dependent upon 
the physical phenomenon represented by the data and the 
desired engineering goals. In this section, the basic con- 
siderations associated with digital simulation are 
discussed . 

Cut-Off Frequency 

Sampling a time history for digital data analysis 
is usually performed at equally spaced intervals of time. 
Care must be taken in determining an appropriate sampling 
interval At. If sampling points are taken too close 
together, they will yield correlated and highly redundant 
data, thus unnecessarily increasing the labor and cost of 
calculations. On the other hand, sampling at points which 
are too far apart will lead to confusion between the low- 
and high-frequency components in the original time history; 
this is so-called aliasing. Generally, at least two 
samples per cycle are required to define a frequency 
component in the data. Hence, the highest frequency which 
can be defined by sampling at the rate of 1/At samples per 
second is l/2At cps. This is called the Nyquist frequency. 

f = l/2At . 


(3-47) 
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Frequencies in the original data above l/2At cps 
will be folded back into the range from 0 to 1/2 At cps, 
and be confused with data in the frequency range below this 
value . 

One way to avoid the aliasing problem is to choose 
At sufficiently small, that is, f large, that it is 
physically unrealistic for data to exist above f . In 
general, it is a good rule (see [17]) to select f to be 
one and one -ha If to two times greater than the maximum 
anticipated frequency. For example, for atmospheric turbu- 
lence simulation, the energy spectrum (as shown in Figure 
2-7, page 19) indicates most of the kinetic energy is in 
the range below 1 cps. Therefore, a meaningful At is about 
0.1 to 0.5 sec . 

The Quantization Error 

Since the magnitude of each data sample must be 
expressed by some fixed number of digits, only a fixed set 
of levels is available for approximating the infinite 
number of levels in the continuous data. For typical 
analog to digital conversion, the quantization error will 
have a uniform probability distribution with a standard 
deviation of approximately 0.29 Ax (see [17], where Ax 
is the quantization increment. This is demonstrated as 


follows : 
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Let P(x) be the quantization error probability 
density function defined by 

P(x) = 1 -0,5 < x < 0.5 

— 0 otherwise . 


The variance of the error 


Since x ?= 0, 


a - 
x 


/ (x-x)^ P(x) dx 


0.5 


x 


= / 


-0.5 


2 , 1 
x dx = Y2 f 


the standard deviation is 

cr = j/ 1/12 - 0,29 scale unit . (3-48) 

X 

In practice, the quantization error is usually 
unimportant relative to other sources of error in the data 
processing. For example, for simulated turbulence, the 
quantization increment Ax is usually chosen as 0.01 m/sec. 
Therefore, the standard deviation of the quantization 
error is approximately equal to 0.0029 m/sec, and is small 
enough to be neglected. However, care must be exercised to 
assure that the range of quantization is small enough to 
assure the accuracy desired. 
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The Bias Error 

All of the power spectral density algorithms con- 
sidered herein use the discrete Fourier transform. The 
finite length of time, T = N*At, influences the spectral 
function. The estimated spectrum <j> is calculated by the 
finite length time history x T (t), which is the product of a 
"boxcar" function b T (t) and the original time history x(t) . 

x T (t) = b T (t) • x ( t ) , (3-49) 

where 

b T (t) = 0 t < T/2 

= 1 -T/2 < t < T/2 

=0 t > T/2 . 


The Fourier transform of b T (t) gives 

°° -i2TTft 

B (f) = /. b (f) e dt 

— QO 


= T ( 


sin (7rfT) * 
irf T 


(3-50) 


this is called a window function. 

Thus, the true spectrum will be different than 
the spectrum <f> estimated by the finite Fourier transfor- 
mation. The two functions are related by the convolution 


integral 
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4> e (f) = s '!> t (f 1 ) B T (f-f 1 )df 1 

^oo 


00 sin ( TT ( f-f , ) T ) 

= T < tf T f-r/y — df x • 


(3-51) 


The difference between <|> e and is defined as the 
bias error. The strongest objection to the boxcar window 
is "leakage through the side lobes," as shown in Figure 3-4. 
The side lobes are those portions of the window between 1/T 
and 2/T, between T/2 and 3/T, and so forth. 



Figure 3-4. The spectral window for boxcar 
weighting function. 


This deficiency in the boxcar window can be over- 
come largely by using a tapered window which tapers the 
ends of the time history. For example, the Hanning window 
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is defined as 

h T (t) — 0 for t < - T/2 

= i I 1 + COS ( ) ] - T/2 < t < T/2 (3-52) 

=0 t > T/2 . 

Then the Fourier transform of h T (t) gives the Hanning 
window function H T (f) where 

H T (f) = \ B t ( f -f Q ) + | B t (f) +i- B t ( f + f o ) , (3-53) 

where 

f = 1/T . 

The Hanning window function is shown in Figure 3-5, 
the side lobes of which are greatly reduced. 




CHAPTER IV 




RESULTS AND DISCUSSION OF SIMULATION 

Chapter IV discusses simulation of turbulence 
employing the various filter techniques described in 
Chapter III. Included in the discussion are simulated 
turbulence time histories, probability density functions, 
power spectrum functions , and coherence functions . Some of 
the simulated results are compared with turbulence 
properties measured in the atmosphere. The atmospheric 
data were measured at NASA Marshall -Space Flight Center, 
Atmospheric Boundary Layer Facility [19] . 

A. Random Number Generation 

A 12-term equation (see Equation (3-6) ) is used for 
generation of the Gaussian random numbers. These numbers 
used as input for the simulations described later are shown 
in Figure 4-1. 

3 

2 

1 
0 

-1 
-2 
-3 

Figure 4-1. Gaussian distributed random signals with zero 
mean and unity standard deviation (2048 data) . 
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By using the computer subroutine GAUSS, as 
described in Chapter III, Section A, and shown in Appendix 
A, Section A-l, ten sets of data, each set having 2048 
digitized points, are generated. The first starting* value, 
IX in subroutine GAUSS, is arbitrarily chosen as 65549. 

The starting value for the successive sets is taken as the 
last output value, IX, from the previous set. The input 
mean value x and the standard deviation a are selected as 
zero and unity, respectively. The properties of each of 
the ten sets of output are shown in Table 4-1. 


Table 4-1. The Variance, a, and Mean Value of Random 
Signal Generated by Equation (3-6) 


Set No. Starting Value Variance cr Mean 


1 



65 

549 

0 . 975 

0.99 

-0.02 

2 

1 

873 

182 

733 

0.962 

0.98 

-0.01 

3 


525 

074 

445 

0.961 

0.98 

-0.02 

4 


250 

707 

981 

0.976 

0.99 

-0.01 

5 

1 

050 

083 

341 

0.983 

0.99 

-0.03 

6 


775 

716 

877 

1.062 

1.03 

0.04 

7 

1 

575 

092 

237 

0.960 

0.98 

0.02 

8 

1 

300 

725 

773 

0.978 

0.99 

0.00 

9 

2 

100 

101 

133 

0.981 

0.99 

0.00 

10 

1 

825 

734 

669 

1.041 

1.02 

0.03 

Average 





0.988 

0.994 

0.00 
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Table 4-1 shows that the generated standard devi- 
ations and mean values are very close to the designed 
values . 

To examine the probability density function of the 
data set, the computer program PDF is used (see Appendix A, 
Section A- 2) . This computer program is based on the 
following equations . 

The probability density function P (x) can be 
defined as (see [17]) 


P(x o ) 


Ax = Prob [ 


/ Ax w ... 
(x o -^-)<x(t) 


, , Ax » , 

( + -7T- ) ] 

O 2 


(4-1) 


where the function Prob [A] represents the probability that 
Statement A is true. More precisely, 


Prob 


P(x o ) 


lim 

Ax->-0 


r / AX, 

[ < x i - — ) 


< x(t) < 


/ L Ax v , 
( x o + - 2 - ) ] 


Ax 


(4-2) 


For digitized data, x n , n = 1,2,...,N, the 
probability that x n falls within the range [ x Q ± x/2 ] can be 
calculated as 


N 


Prob [ ( x Q - ) 


A X 

^ / , Ax , o 

X < ( X^ + -77- ) ] = —rr~ 

n - o 2 N 


(4-3) 


where N„ is the number of 
x o 

the range [ x q ± Ax/2 ] . 


digitized values which lie in 
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The probability density function defined in 
Equation (4-2) can thus be written as 


P(x ) = lim . (4-4) 

Ax+O “ flx 

Rather than using the limiting process for digital 
simulation, a small, finite value of Ax is used which gives 


N 

- x o 

P<X o ) - N Ax 


(4-5) 


Figure 4-2 shows the calculated P(x q ) for Data 
Set 1, using Ax = 0.3, compared with the theoretical 
Gaussian distribution defined in Equation (3-4) . Figure 
4-3 shows the same P(x q ) computed by using a Ax value of 
Ax = 0.2. Comparison of the two figures shows that the 
variation in Ax has little influence on the agreement of 
the P (x q ) function with the theoretical distribution to be 
simulated. Figure 4-4 shows the influence of using 
Ax = 0.2 and of increasing the sample size, N. Using the 

A 

total 20480 datum points, the estimate of P (x q ) agrees con- 
siderably better with the actual Gaussian distribution. 

As mentioned previously, a Gaussian distributed 
white noise input is required for the proposed simulations. 
The power spectrum for white noise is constant for all 
frequencies. The power spectrum for the generated random 
signal input was computed to confirm that it was, indeed, 


constant . 


Figure 4-2. The estimated probability density function 
of the 2048 random numbers, shown in 
Figure 4-1, compared with the theoretical 
Gaussian distribution (Ax = 0.3). 
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A raw estimate of the power spectrum is defined as 
(see 117] ) , 


((.(f) = | |X(f,t) | 2 , 


(4-6) 


where T = N*At, and X(f,T) is the Fourier transform of x(t) 
for discrete frequency values shown in Equation (3-15) , 

X ( f m' T > 

X m = ST m = 0,1,2,...,N -1 . (4-7) 


Hence , the discrete power spectrum becomes 


<i> 


m 



(4-8) 


To transform the time history x^ to the frequency 
domain X , the fast Fourier transform computer subroutine 
FFT (see Appendix A, Section A-3) is used. The estimated 
discrete power spectrum is then calculated from Equation 
(4-8) . The computer subroutine for this purpose is called 
SPEC (see Appendix A, Section A-4) . The estimated spectrum 
of Data Set 1 as generated by this subroutine is shown in 
Figure 4-5 (a) . 

Since the estimated spectrum is for bandwidth- 
limited input data, estimates at a frequency spacing of 1/T 
will be essentially uncorrelated. Hence, smoothing or 
averaging techniques are required to obtain a representa- 
tive spectrum function. One smoothing technique consists 



ii . I, Li i>» .?■ X Jl . \ AjU . .lli W* u » J ik 1 ' ■* u» l*_k?jl4^Al j_ ^jbj-i-±_jL 


(a) Estimated raw spectrum 


1/ (20At) 


(b) Ten-point frequency smoothed spectrum 


(c) Segment-averaged spectrum for ten data sets 

Figure 4-5. The estimated raw and smoothed spectra of 
Data Set 1 (Figure 4-1, page 54) 

(At = 0.1 sec) . 
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of averaging n neighboring frequency components of the raw 
spectrum estimates (see [17]), 

= I 1 + Vl + • ' ' + W-l 1 • (4 ' 9) 

Figure 4-5 (b) shows the frequency smoothed value <J>^ of the 
<J> given in Figure 4-5 (a) . 

A second technique of smoothing is averaging the 
results from n separate time slices, segment averaging, 
where each slice is of length T such that the original 
length is equal to n«T. The averaged spectrum is given by 
(see [17]) 


K = w f 4> m i + <f> 

m n m, 1 


m, 2 


+ 


+ 


<J> ] , 

r m,n ' 


(4-10) 


where the second subscript indicates the slice number. 

Figure 4-5 (c) shows the segment-averaged values of 

d)", as well as the average of the <b values calculated 
m m 

separately for each of the ten sets of data listed in 
Table 4-1. The computer subroutine for smoothing is 
called SMOOTH and is listed in Appendix A, Section A-5. 

B. Simulated Turbulence 

As discussed in Chapter II, various spectrum models 
can be simulated by different filter functions. A general 
form of H(f), which represents the Dryden filter and is a 
good approximation of the Kaimal and von Karman spectrum 
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filter, is given by 


H (f ) 


c + jdf 
[ a + jf] * 


(4-11) 


where the constants a, c, and d are different for each 
model, as listed in Table 4-2. 

Equation (4-11) can be written as 


H(f) 


c + jdf 

[ ( a 2 - f 2 )+ 2 jaf 

D [ c ( a 2 - f 2 ) + 2adf 2 + j [ ( a 2 - f 2 ) d - 2ac ] f ] , 


where 


(4-12) 


D ~ 


("a 2 -f 2 ) 2 + 4a 2 f 2 ‘ 


The real and imaginary parts of H(f) function are then 
given by the two functions 


HR (f ) - Die ( a 2 ^ f 2 ) + 2adf 2 ] 


HI (f ) = D [ ( a 2 - f 2 ) d - 2ac ] f 


(4-13) 


A computer subroutine FILTER (see Appendix A, Section A-6) 
is available for the calculation of the real and imaginary 
parts of the filter function. This subroutine also calcu- 
lates the output Y (f ) from 
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Table 4-2. Coefficients a, c, and d for 
Different Spectrum Model 


a 


b 


c 


Longitudinal component v g V 
of Dryden spectrum —(j) —(j) 
(Equation (2-11)) 




1/2 


Lateral component 3/2 

of Dryden spectrum 1 /A 3(3 /• V 

(Equation (2-11)) 2irt/r ^2 2A 


tt l A j 


1/2 


Approximate v 3/2 1/2 

Kaimal spectrum 0.155(V) 0.0432 a(^) 0.279 C(V) 

(Equation (2-20) ) 


Approximate v 3/2 v 1/2 

von Karman spectrum 0.286(-^-) 0.114 0.398 a(-^) 

(Equation (2-23)) 
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YR(f) = HR (f ) XR (f ) - HI (f ) XI (f ) 
YI(f) = HR (f ) XI (f ) + HI(f)XR(f) 


(4-14) 


where XR(f) and XI (f) are the real and imaginary parts of 
the Fourier transform of the input x(t) . 

The frequency function Y(f) is transferred to the 
time domain, y(t) , with the inverse FFT method. 

The difference equation based on the 
Z-transf ormation technique can be derived for the filter of 
the spectrum given by Equation (4-11) , as illustrated for 
the simple case in Chapter II, Section B. Thus, the 
discrete output Y n+ ^ is related to y f Y n _^/ x n / ant ^ x n -l 
by: 


y n+l = c l’ y n + c 2 * y n-l + d l " x n + d 2‘ x n-l ' 

where 

0 -^ = 2 exp [ -2aTT At ] 
c 2 = - exp t -4ar At ] 

d 2 - d 1^2 < S= S^> it] C l /2 * 


(4-15) 


The constants a, c, and d are those listed in Table 4—2, 


= (2tt) 


1/2 


and the constant b 



66 

A computer program employing the Z -transformation 
technique called DZT is given in Appendix A, Section A-7 . 

The FFT/FILTER method and the DZT method are now used to 
simulate turbulence time histories having the different 
spectra described above. 

Figures 4-6 and 4-7 show the simulated time 
histories for the four spectra as defined by the constants 
listed in Table 4-2, using both the FFT and DZT methods, 
respectively. All data shown in Figures 4-6 and 4-7 are 
normalized with o to facilitate comparison. Smaller values 
of cr DZT therefore give larger peaks, as shown in Figure 4-7. 

Although the expected values of x and a are zero 
'and unity, respectively, the simulated mean value x and 
standard deviation a are affected differently by the FFT 
and DZT methods. Tables 4-3 and 4-4 show values of x and a 
determined statistically from the simulated turbulence time 
histories using Table 4-1, page 55, as the input. Table 
4-3 shows that the value of x is approximately two times 
the input value for the FFT simulation, but almost exactly 
the input for the DZT method. The overall average value, 
however, is still in the acceptable range. 

Table 4-4 indicates that the input variance is con- 
siderably altered by both techniques. This effect is 
strongly dependent on the time increment At. 

To generate the correct variance, one must consider 
the difference between a continuous system before and after 
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(b) Lateral component of Dryden spectrum model 
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(c) Approximated Kaimal spectrum model 
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(d) Approximated von Karman spectrum model 

Figure 4-6. The normalized time histories of simulated 

turbulence by the FFT method using Data Set 1 
as input (At = 0.5, N = 2048) . 
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Table 4-3. The Mean Value of the Simulated Turbulence 



Input 

DZT 

FFT 


-0.02 

-0.02 

-0.06 


-0.01 

-0.01 

-0.02 


-0.02 

-0.02 

-0.05 


-0.01 

-0.01 

-0.02 


-0.03 

-0.04 

-0.09 


0.04 

0.04 

0.09 


0.02 

0.02 

0.06 


0.00 

0.00 

0.00 


0.00 

0.00 

0.01 


0.03 

0.03 

0.07 

Average 

0.00 

-0.01 

-0.01 



Table 4-4. 

The Standard Deviation 
Simulated Turbulence 

a of the 


At 

0.1 

0.2 

0.5 

1.0 

2.0 

Input 

0.99 

0.99 

0.99 

0.99 

0.99 

DZT 

0.18 

0.25 

0.40 

0.56 

0.79 

FFT 

0.44 

0.63 

0.99 

1.36 

1.86 
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it is passed through a sampling and holding device. The 
sample and hold function x(t) is related to the continuous 
function x(t) by 


x(t) = £ x(nAt) f s ( t-nAt) - s ( t - (n+l)At ) ] , (4-16) 

n =0 


where s(t) is the step function defined 


s(t) =1 t > 0 


= 0 t < 0 , 


Taking the Fourier transform of Equation (4-16) 


X(u>) = 


l * n i 

n =0 n 


-jnAtw _ -j(n+l)Atw 


= l X e"3 nAt “ [ 

L r, n 


1 - e _jAt “ 


] 


n =0 




! __ 

= X [ q- ] , 

m 3 w 


(4-17) 


where X is the discrete Fourier component of x(t) (see 
m 

Equation (3-16)). Equation (4-17) therefore becomes 


X (w) = X (w) • s (w) 


(4-18) 


where 


1 r 1 - exp ( - j Atcx)) 1 
At, [ 3^ J 


s (w) 
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Let the continuous function input have a constant 

. 2 

spectrum $ (f) and a variance of a . Then the spectrum of 

X X 

the sample and hold output is 


$ x (£) = • |s(f) I 2 . 


(4-19) 


Since the constant spectrum <J> x (f) can be replaced by 


9 l/2At 

= / 4>„(f) df 

x o x 


4> x (f) = 2At o* , 


it follows that 


$ (f) = 2 At a 2 i »in(Tf|At). ] . 
Y x x 7rf At 


(4-20) 


Using this sample and hold data as the input of the 
turbulence simulation system, the simulated output variance 
is 


a 2 = / $ (f) I H (f ) | 2 df , 
y o 


(4-21) 


and since the continuous system output <f> is defined by 
2 

<f> =* | H (f ) | , Equation (4-21) becomes 


a2 r O 2 A4 _ r Sln(TTfAt) -J , -.r; 

a = / 2a At [ — — *bnr J 4>(f) df , 

y x irfAt y 


(4-22) 


o 
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Assuming { sin(TrfAt) / irfAt ] to be approximately unity for 
the lower frequency (see Table 4-5) where the atmospheric 
turbulence has the most significant energy density, 
Equation (4-22) reduces to 


2At S <p (f)df 
o y 


= 2At a' 


If a x is equal to unity, then 


r<2 2 

a = 2At , 

y y 


or if a is selected equal to 

X 


2 _ 1 
°x 2At r 

then 


Table 4-5. Values of I sin(7rfAt) / nf At J for 
Low Frequency f (At = 0.1) 


0.001 0.01 0.1 


1.0 2.0 


r sin (TTf At ) , 
L TTf At J 


0.999 0.999 0.999 0.995 0.983 0.935 0.858 
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• • 2 
Hence, by adjusting the input variance cr through 

the relationship 1/2 At, the discrete system will have the 

desired output variance a . 

Considering now the difference equation of the 

Z-transf ormation, the output variable in terms of the input 

variables is 


n 


y n - 


I 


i=l 


x i 


(4-27) 


Since each value of x ^ (or y^) is independent, the variance 
can be given by 


a 2 

y 



lim 

n+°° 


V 2 

Z a i 

i=0 1 


(4-28) 


2 

To find the ou , one can write the initial terms of 
the difference equation for the first few inputs and then 
deduce the functional relationship between and 
The sum of the squares of these terms can then be split 
into geometric series and derivatives of geometric series, 
which are all summable. For example, the result for 
Equation (3-46) is (see Appendix B) 


^2 

a 

y 


= a . 


x 


2a 2 

[ y 


( 1 - e -aAt ) 


Tra 


1 - e 


-2aAt 


(4-29) 


where 


V 

a “ X • 



For small At, Equation (4-29) reduces to 


s 2 

JL 

a 2 

X 


M a 2 

ir y 


With a = 1 the result is 


-2 At 2 

a = — a 
y 7T y 


(4-30) 


(4-31) 


Alternately, if the input is selected as 


2 TT 

G x “ At ' 


(4-32) 



2 

be equal to the desired output variance a^.. 

Figure 4-8 compares Equations (4-24) and (4-31) 


with the calculated standard deviations given in Table 4-4, 
page 69. In both cases good agreement with the theory in 
the small At range is observed. This is expected, since 
Equations (4-24) and (4-31) also assume small At. Hence, 
the variance can be simulated properly by applying 
Equations (4-25) and (4-32) , respectively, to the input of 
the system. The computer program for generating the input 
array x^ with the adjusted variance is called INPUT, and is 
listed in Appendix A, Section A-8. 

To compare the simulated turbulence with measured 
atmospheric turbulence, the value of the reference velocity 
V, the turbulence integral length scale A, and the 
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turbulence intensity a must be specified. The reference 
velocity is usually selected as the mean wind speed at the 
given height of the simulation. The A and a values are 
functions of height, z, atmospheric stability conditions, 

( z/L) , and surface roughness, z q . The parameter L is the 
Monin Obukhov stability length. Equations for predicting 
the turbulence intensity values are given in 17] as: 


o, = 0.52 U/ I An ( + 1 > + Mz/L> 3 

o 

o -0.8 

a 2 = a 3 I 0.583 + 1.39x10 z ] 

o -0.4 

a l ^ a 3 I °* 177 + 2.74x10" J z ] 


(4-33) 


The turbulence length scale for the Dryden spectrum is 
given in 116] as; 


A 


1 



A 


3 


44.21(3.28 z) 1/2 

533 

z 


z < 533m 


z > 533m 
z < 533m 


533 


z > 533m . 


(4-34) 


Figure 4-9 shows the comparison of the longitudinal 
velocity component time histories of simulated and measured 
atmospheric turbulence. The mean wind speed, U, and 
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turbulence intensity, a , of the measured data are 7.52 m/s 
and 1.78 m/s, respectively. Based on this value and a 
height of 24 m for neutral stable conditions, the calcu- 
lated value of from Equation (4-33) is 1.26 m/s, which 
is the standard deviation used to simulate turbulence. 
Hence, the time histories shown in Figure 4-9 are 
normalized to facilitate comparison. From visual obser- 
vation of the time histories, it appears that the measured 
data contain more high-frequency components than the simu- 
lated results. This is expected because the Dryden form of 
the spectrum results in a filter which generates less 
power at higher frequencies (see Figure 2-2, page 9). 

The probability density functions (PDF) of the 
measured and simulated turbulence are given in Figure 4-10. 
All of the simulated turbulence is very nearly Gaussian 
distributed; however, the measured atmospheric turbulence 
has non-Gaussian characteristics, as discussed in Chapter 
II, Section C. Hence, the multi-filter system shown in 
Figure 2-8, page 20, and in Table 2-1, page 26, has been 
programmed (Appendix A, Section A- 9) to provide better 
simulation of the non-Gaussian characteristics of atmos- 
pheric turbulence . 

The calculated PDF’s of this model for different 
values of r, the adjusting coefficient, are shown in 
Figure 4-11. A proper selected value of r is observed to 
provide a correct PDF for the simulated turbulence. 
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J 

0 




(a) Measured atmospheric (b) Simulated Dryden spectrum 

turbulence turbulence 



Simulated turbulence 
of approximate 
Kaimal model 


(d) Simulated turbulence of 
approximate 
von Karman model 


Figure 4-10. The estimated PDF of the data of the 
different spectra techniques 
(Ax = 0.2) . 


(c) 
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Figures 4-12 and 4-13 show the time histories and PDF's of 
the simulated non-Gaussian turbulence and of the measured 
atmospheric turbulence, respectively. 

The PDF of the non-Gaussian model gives a peak shape 
similar to the measured data. This model also increases 
the number of sharp peaks relative to the Gaussian model; 
however, the number of overall high-frequency components is 
still much less than that for the measured atmospheric 
turbulence. This is to be expected because the non-Gaussian 
model still utilizes the Dryden spectrum filter. 

In order to examine how well the simulation repro- 
duces the input spectra functions, the discrete estimates 
were calculated from Equation (4-8) . The segment-averaging 
technique. Equation (4-10) , is applied to smooth the 
spectra . 

The spectra estimates of turbulence are shown in 
Figures 4-14 and 4-15 for the FFT and DZT methods, 
respectively. Each plot is smoothed using 20 sets of data. 
The theoretical curves, based on Table 4-2, page 64, are 
also plotted for each case. Comparing these two figures 
shows that the FFT method provides a better spectrum simu- 
lation than the DZT method. 

The spectrum of measured atmospheric turbulence was 
also computed and is compared with the von Karman spectrum 
(Equation 2-10)) and the Dryden spectrum (Equation (2-11)) 
in Figure 4-16. In this particular case, the measured 
spectrum agrees best with the von Karman spectrum. 



• T J .Ml 


(a) Measured atmospheric turbulence (z 


100 sec. 


ril 


LTVi*] 


(b) Simulated non^-Gaussian Dryden spectrum model 
Figure 4-12. The normalized time histories (r = 1.5 


(a) Measured atmospheric 
turbulence (z — 24 m) 


(b) Simulated non-Gaussian 
Dryden spectrum model 
(r =* 1.5) 


Figure 4-13. The estimated PDF of the non Gaussian 
turbulence . 



83 


< 


Sf 

<y 



(by FFT method, At = 0.5, N = 2048 
Df ten sets) . 




(c) Approximate Raima 1 (d) Approximate von Karman 



-14. (continued) 





(by DZT method, At = 0,5, N = 2048 
of ten sets) . 






Figure 4-15. (continued) 





Figure 4-16. 


Comparison of two theoretical curves with the 
calculated spectrum of measured atmospheric 
turbulence (z = 24 m f At = 0.5, N = 2048) . 
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It is well known that the von Karraan spectrum more 
accurately describes atmospheric turbulence than does the 
Dryden spectrum (see Figure 2-2, page 9). For stable 

atmospheric turbulent flow conditions the Kaimal spectrum 
is required , 

Because of the irrational nature of the von Karman 
and Kaimal spectrum functions, no real physical filter 
function can be obtained for simulation purposes. However, 
an exact solution of the filter function for these two 
spectra can be obtained. For the von Karman spectrum, the 
filter function is given as: 


- : 7 . . 5/6 

(a + ju) 7 


(4-35) 


H 2 (w) = 


b(juj + c 2 ) 

(a + ' 


where 


V 


a = (0.74682)^ 


c, = (0.62558) [ 


b = 


a 2 V 1//3 

(0.72236 [ - l7 y -j 


A' 


a l V 


1/3 


A 


1/3 


] , c„ = (0.4573) 


V 

A * 


For the Kaimal spectrum the filter section is: 


c 

a + j (f ) 


576 


r 


H(f) 


(4-36) 
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where 

V tt 1/3 

a = (0.17241) [~] , c = (0.34482) [ £ ] 

Therefore, by using Equations (4-35) and (4-36) in con- 
junction with the EFT method as described in Equations 
(4-13) and (4-14) , the exact von Kantian or Kaimal turbu- 
lence spectra can be simulated. The computer program for 
simulation of these two spectrum functions is included in 
the subroutine FILTER (Appendix A, Section A-6) . 

In Figure 4-17 the time histories utilizing these 
more exact spectrum functions are compared with measured 
atmospheric turbulence. It can be seen from this figure 
that simulation of the high-frequency components is 
improved greatly over the Dryden and approximate filter 
simulations shown in Figure 4-9, page 77. The spectrum 
functions computed by Equation (4-8) for these two cases 
are shown in Figure 4-18, and the results agree very well 
with the theoretical curves. 

Two disadvantages for these mathematical models 
are: first, due to the irrational form of the filter 

function, the DZT method cannot be used to establish the 
simple difference equation; and second, the correlation 
functions for both the von Karman and Kaimal spectra are 
much more complicated than the exponential function for the 
Dryden spectrum (see Equation (2-26)). For example, the 
correlation function for the von Karman spectrum is in the 
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(a) von Karman spectral (fr) Kaimal spectral 

model model 

Figure 18. Comparison of theoretical with simulated 
spectrum ( At = 0 . 5 , N = 2048 , segment- 
average of ten sets) . 
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form of modified Bessel function of the one-third order [5] . 
Thus, the von Karman and Kaimal spectra cannot be utilized 
for the multi-filter system to simulate the non-Gaussian 
character of atmospheric turbulence. 

C. Simulation with Interlevel Coherence 

In this section, simulated horizontal components of 
wind speed which include vertical coherence are discussed. 
Measured wind data recorded at the eight-tower Atmospheric 
Boundary Layer Facility, Atmospheric Sciences Division, 

NASA Marshall Space Flight Center, are also computed and 
presented for comparison. 

The output transform, discussed in Chapter II and 
in reference [3] , is given by 

M 

Y(K,z) = A H(K,z) l D m (K,z)X m (K) , (4-37) 

m=-M 

where the D m (K,z) was defined by Equation (2-38) as: 

V K ' Z) = c m exp I jKzmir / £ o ] , = (K Az) max . 


Substituting D^(K,z) into Equation (4—37) gives 


Y(K,z) 


M j 9 

H(K,z) [ l ( e m X m (K) 
m=0 


-j 0 


m 


X 


-m 


(K ) ) B 1 


m 


+ e 


(4-38) 


where 
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c A 

B 0 = -V" ' B m = c m A for m * 0 

0 = xnKz7T / e 

m o 

Equation (4-38) can be reduced to 

Y(K,z) = H(K,z) X'(K,z) , (4-39) 

where the 

M j 0 -j0 

x' (K,Z) = l ( e “ X (K) + e m X (K) ) B . 
m=0 

Thus, the transform of the output Y(K,z) can be calculated 
in the same manner as described in Equations (4-13) and 
(4-14) . The computer program for calculation of these 
equations is listed in Appendix A, Section A-10 . 

Simulated results are obtained by using 11 levels 
of random signal inputs (M equal to 5 in Equation (4-37) ) 
and the Dryden spectrum filter function for H(K,z) . 

Figures 4-19 and 4-20 compare the time histories and auto- 
spectra of simulated coherent turbulence with the measured 
values for the 24 m and 12 m levels. Both figures show 
simularly shaped curves of equal magnitude. 

As also mentioned in Chapter II, the atmospheric 
turbulence near the ground shows a coherence which behaves 


as 





0.01 o.l 1.0 10.0 


001 0.1 l.o 10-0 


(a) Simulated interlevel coherence turbulence at 
z equal to 24 m and 12 m, respectively 


0.01 0-1 10 10.0 


0-01 0-1 1-0 10.0 
f 


(b) Measured atmospheric turbulence at z equal to 
24 m and 12 m, respectively (At = 0.5) 

Figure 4-20. The spectrum estimated of level 
coherence turbulence (N — 2048) . 
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Y 0 (n> - exp [-an] , 

where 


(4-40) 


Az f 

n = “v"" • 


This expression for Y 0 (n) was originally proposed by 
Davenport [13] as a simple scaling law to describe the 
coherence of the vertical interlevel winds. He found a to 
be equal to 15.4 (note some reports use 7.7 for the square 
root coherence function) . Stegen and Thorpe [12] show this 
is approximately correct. Work by Pielke and Panofsky [21] 
and Brook [11], however, provides values for the coeffi- 
cient, a f which are dependent on stable conditions, as 
shown in Table 4-6. 


Table 4-6. The Estimates of the Coefficient, a, of 
Different Stable Condition [11] 


Component 


n Range 


Stable 


Neutral 

Stable 


Unstable 


Longitudinal 


Lateral 


all - n 
r\ < 0.12 


all - n 

T) < 0.12 


16.8 

13.4 


12.5 

14.4 


17.0 

13.2 


12.8 

13.4 


16.8 

16.8 


10.4 

11.6 
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The coherence function is usually calculated from 
spectra as follows (note: here the x and y simply indicate 

two different time histories) : 


I <j> (r \ ) P 

Y xy (n) = ' 

where the cross -spectrum is defined as: 


(4-41) 


<J> xy (n) = ~ [X(n) Y*(n) ] . 


(4-42) 


By using this equation and Equation (4-8) for the auto- 
spectra, the coherence function for a multi-input system 
as shown in Equation (4-37) can be calculated as: 


Y (n) 


2 M 2 : 

1 c 0 + 2 2 c m COS ( mnl,/ Vx 3 ] 

m=U 

~ M 

t c 0 + 2 I' 4 3 

m=0 


(4-43) 


where 


max 


A. z f 

max max 

V 


The value of f is taken as the cut-off frequency 
max 

(Nyquist frequency) for the discrete transformation. In 
order to reproduce accurately the coherence function y q 
given in Equation (4-40) , the value of ^ z max mus ^ 
selected appropriately. Figure 4-21 shows the exponential 
curve with a = 17 compared with Equation (4—43) for 
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different values of Az max - This figure indicates that 
small values of Az max , approximately equal to 2 or 3, 
represent reasonable values to describe the exponential 
curve. Note that Equation (4-43) produces a second peak 
after one period at 


n / n 


max 


2.0 


(4-44) 


therefore, in programming the coherence model, the 
undesired cyclic peaks occurring at high frequencies must 
be truncated, A critical frequence n c is defined by 


n 


c 


0 


max ’ 


(4-45) 


This critical frequency defines the upper limit of the 
coherence function in Equation (4-43); hence, 

y(n) = y(n_) for all n > n . (4-46) 

G G 

Figure 4-22 compares the calculated coherence 
functions of simulated and measured turbulence with the 
exponential curve given by Equation (4-40) . This figure 
shows that the model behaves very much like the atmospheric 
turbulence and, at the same time, indicates that the calcu- 
lated coherence functions are described appropriately by 
the exponential curve . 
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The coherence of the simulated time histories was 

computed using statistical techniques. Considerable care 

must be taken in evaluating the coherence function because 

the statistical estimate has a rather complex sampling 

distribution. The coherence estimate must be smoothed by 

either ensemble or frequency averaging in order to suppress 

the random error to acceptable levels. Otherwise, even for 

the case of totally incoherent data, the estimated 

coherence will be unity for the entire frequency range. It 

should be noted, in particular, that the smoothing must be 

carried out separately for the real and imaginary terms of 

the cross-spectrum shown in Equation (4-42) . Therefore, 

<J> (n) can be represented by 

xy 

^xy (n) = ^xy (n) + ^xy (n) ' (4-47) 


where the <j>^ and (J>^ are the real and imaginary parts of 


xy r 


called the co-spectrum and the quad-spectrum. 


respectively. To smooth the cross-spectrum, the averaging 
technique must be applied to each term individually. 
Therefore , 



(n) = <|> 1 

xy 


<n) + j^y (n) 


(4-48) 


where the bar indicates the averaged value. Then the 
coherence Y X y( r l) can be calculated by 
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y (n) 

? xy 


l , f > xy < T i) l 2 
(n) ?__(n) 

Y 


( 4 - 49 ) 


It is re-emphasized that the average must be taken 
before the absolute function of the cross-spectrum is 
computed. This is because the value of the spectrum 
function is small relative to the high variation created by 
the discrete transformation. Therefore, the magnitude of 
the standard deviation of the variate is much greater than 
the magnitude of the mean. Thus, smoothing the absolute 
value of the spectrum is, in effect, just averaging the 
variance. The variance is always positive; therefore, the 
average is very large. By averaging <j>^. and <j>^. 
separately f the positive and negative values reduce the 
magnitude of the average. The computer program for calcu- 
lation of the coherence function and the smoothing 
techniques is listed in Appendix A, Section A-ll. 



CHAPTER V 


APPLICATION AND CONCLUSION 

The purpose of this chapter is to describe how the 
different turbulence simulation techniques influence 
computer-simulated landings of aircraft having the charac- 
teristics of a DC-8. A three-degree-of-f reedom dynamic 
model of an aircraft during approach is utilized with 
different models of simulated turbulence imposed as atmos- 
pheric gusts . The development of the model of the airplane 
dynamics and a description of the computer program are dis- 
cussed in detail by Frost and Reddy [22] and Frost and 
Crosby [23] . 

The three-degree-of-freedom airplane computer model 
uses Runge-Kutta integration to solve the system of time- 
dependent, second-order ordinary differential equations 
(see Appendix C) . The flight simulation program incorpo- 
rates an automatic control system in which the basic 
control inputs are thrust and pitch angle and the control 
variables are speed and flight path height. The model 
selector automatically selects the proper control modes in 
sequence according to the predetermined flight path. The 
four modes are hold, capture, track, and flare, as shown in 
Figure 5-1. 


103 


104 


altitude 
hold mode 



Figure 5-1. Automatic landing control modes. 

A. Method of Inputting Wind 

For input to the dynamic equations , the actual wind 
speed U at location (x,z) and time (t) is calculated by 

U i (x,z,t) = U ± ( z ) + u ± (x,z,t) , (5-1) 

where the mean wind speed, U, upon which the turbulence is 
superimposed, is described by a logarithmic function of 
altitude z. Thus, 

u * z + z _ 

D,(z) = In I — — ° ] , U 2 = 0 , (5-2) 

o o 

where 

u* is the friction velocity 

K is the von Karman constant 
o 

z Q is the surface roughness length . 
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To utilize the simulated turbulence described in 
the preceding chapters for analysis of landing aircraft, 
u i (x,z,t) is replaced by the discrete signals. Inputting 
the discrete data into the computer program requires some 
special attention: 

1 . Simulation calculates N values for a given 
height z. 

2. The airplane is descending, so z is changing 
with each time step. 

3. A set of signals from a given transformation is 
generated for a period of 12.8 seconds, which 
is selected by using N equal to 128 (2 ) and 

At equal to 0.1 second. 

It is assumed that the turbulence characteristics 
encountered by the airplane during that time period are 
uniform. However, as a second control, the program con- 
tinuously checks altitude. If during the 12.8-second time 
period the aircraft has a sudden altitude change (in this 
case, 20 percent of the height z) , the program auto- 
matically generates a new set of turbulence signals at the 
new height. 

A second consideration is that the simulated turbu- 
lence is the discrete time history of wind fluctuations at 
a fixed point, and Taylor's hypothesis must be applied to 
establish a space and time relationship for calculating the 
effects due to location change, as shown in Figure 5-2. A 
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(a) A fixed point P, after 
7At time steps, will 
encounter the eighth 
turbulence signal 


(b) A moving point P, after 
7At time steps, will 
encounter the thirteenth 
turbulence signal 


Figure 5-2. Effects of aircraft motion relative to atmospheric motion 
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time scale, t^. defined in Equation ( 5 ^ 3 ) , is therefore used 
to count the turbulence signals 

t L = t + ~ = n At , (5-3) 

where 

t is the time scale for counting turbulence signals 

t is the real time measured from the airplane entering 
the wind field 

d is the distance moved by the airplane in time t 
V is the airplane relative velocity . 
th 

Hence, the n one of N-generated data is chosen as the 
turbulence fluctuation u^ at the time t. The lift and drag 
(see Appendix C) on the airplane are therefore computed 
from the wind UN , defined as 

U.(t) = U. + u. (n At) . (5-4) 

l ii 

The wind also enters into the equation (as shown in 
Appendix C) by 


dU. 

l 

dt 




9x 
9 1 


+ 


30. 
3 z 


9 z 

9 1 * 


(5-5) 


Hence, each term becomes two parts. For example, the 
9tN/9x becomes 


9U i 91K 3u i 
9x ~~ 3x + 9x 


(5-6) 


I 
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Because the Taylor hypothesis is used, the gradient 
S^/Dx of Equation (5^6) is represented by 

3u a 

= lu. ( (n+1) At ) - u ± (n At) ] / (V At) . (5-7) 

The computer program which provides the wind field 
data and the gradients for the aerodynamic equations is 
listed in Appendix A, Section A-12. Since the vertical 
motion of the airplane is relatively small compared with 
the horizontal motion, the gradient terms in the Z - 
direction can therefore be neglected. 

B. Simulated Landings 

The simulated landing begins at an altitude of 300 
meters and a distance of 7 km from the runway in the 
altitude hold mode. Although this height is somewhat low 
for an actual landing case, it has been used to save 
computer calculation time. The glide slope angle is 
2.7 deg and the initial trimmed airspeed is 70 m/s. 

The major aircraft landing error is the touch-down 
position shifting, which is mainly affected by mean wind 
shear and turbulence. The mean wind shear is dependent 
upon the parameters u* and z q , as described in Equation 
(5-2) . For a no-turbulence case, the simulated results of 
landing positions are plotted against different values of 
u* and z q in Figures 5-3 and 5-4, respectively. A head 
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wind decreasing vertically toward the ground causes the 
plane to land short of the runway, whereas a tail wind 
decreasing vertically toward the ground causes the plane to 
overshoot the landing point. The scatter in results as 
seen in Figures 5 ->3 and 5-4 is probably due to the over- 
response of the automatic control system. 

In the landing simulation in which turbulence is 
the main cause of landing position error, different turbu- 
lence models have been applied. The statistical result for 
each model is based on a sample of 20 simulated landings. 
For the simple turbulence model with a Dryden spectrum, the 
statistical results of landing position are shown in 
Table 5-1, 

The turbulence intensity, a, strongly influences 
scatter in the landing position about the mean; i.e., 
landing position error. In the discrete transformation the 
output standard deviation, a, can easily be adjusted by 
using Equation (4-23) or Equation (4-30) to give a constant 
ratio of the original calculated turbulence intensity. It 
is not realistic of changing turbulence intensity without 
changing the value of u* . The assumption has been made for 
only showing the influence of different values of o. 

The error in the average landing position with 
respect to the no-turbulence condition, defined as the 
average error, and the standard deviation for a plane 
landing in Gaussian and non-Gaussian turbulence with a 


Ill 


If 


Table 5-1 . Twenty Landings for Simple Turbulence Model 
using a Dryden Spectrum with u* = 0.5 and 
z = 0.1 (FFT Method) 


Distance from Average 

Run No. Landing Position Touchdown Point 


1 

6919.07 

30.71 

2 

6883.52 

-66.26 

3 

6950.17 

0.39 

4 

7000.64 

50.86 

5 

7001.87 

52.09 

6 

6969.55 

19.77 

7 

6988.47 

38.69 

8 

6992.39 

42.61 

9 

6907.23 

-42.55 

10 

6948.49 

-1.29 

11 

6956.65 

6.87 

12 

7070.36 

120.58 

13 

6970.75 

20.97 

14 

6945.24 

-4.54 

15 

6900.77 

-49.01 

16 

6922.49 

-27.29 

17 

6910.08 

-39.70 

18 

6891.43 

-58.35 

19 

6893.70 

-56.08 

20 

6972.73 

22.95 

Average 

6949.78 S.D. 

46.37 

Note: The results for non-Gaussian, Dryden 

spectrum, and for the Gaussian, von Karman spectrum 
turbulence landing simulation are listed in Appendix D. 
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Dryden spectrum, and in Gaussian turbulence with a von 
Karman spectrum, are shown in Figures 5-5 and 5-6. Both 
figures show the direct relation of turbulence intensity to 
deviation in landing position. 

Another important factor relative to landing per- 
formance is the touchdown sink rate. From Figures 5-7 and 
5-8 one can conclude that the turbulence intensity has 
little influence on the average touchdown sink rate for the 
automatic control system used in this analysis. However, 
the variability of sink rate, i.e., standard deviation, 
does increase with increasing turbulence intensity. This 
occurs because the touchdown sink rate is very sensitive to 
the turbulence fluctuation after the flare mode is started. 
Near the ground U becomes small , and if an increased head 
wind due to turbulence is experienced followed by a large 
reduction, or even a tail wind near the runway, a sudden 
loss of lift occurs and causes a very hard landing. Also, 
as noted in Figure 5-8, when o is over 3 to 4 m/s ec, the 
standard deviation of sink rate approaches 30 percent of 
the mean sink rate. Thus, the chance of a hard landing 
becomes high. 

The influence of different turbulence models for 
this landing simulation can also be seen in Figures 5-5 
through 5-8. The average values of position error and sink 
rate are affected little by using different turbulence 
models, whereas the standard deviation in touchdown 
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position is influenced significantly by the different 
turbulence models. 

For the same spectrum function, a non-Gaussian 
turbulence model produces the higher degree of touchdown 
variability. It should also be noted that , when different 
spectrum functions are used, the von Karman model induces 
a smaller standard deviation in position error than the 
Dryden model. These flight simulation results indicate 
that for a DC-8-type airplane, the turbulence energy con- 
tained by the higher frequency fluctuation affects the 
landing position little. The sharp peaks of the non- 
Gaussian model represent the most dangerous events which 
can cause the airplane to lose control during the approach. 
Hence, it is important to simulate properly the non- 
Gaussian nature of the atmosphere. 

The DZT method for turbulence generation and the 
interlevel coherence turbulence model were also used to 
carry out flight simulation. As expected, the statistical 
results are essentially the same for the DZT method as they 
are for the FFT method. Also, the interlevel coherence 
model was found to give approximately the same results as 
the one^-point spectrum model. The reason for this is that 
the airplane, when flying along the 2.7-deg glide slope, 
travels relatively longer distances in the horizontal 
direction than in the vertical direction. The coherence 
between vertical levels has little effect on the two-point 


1 1 6 


separation over relatively large horizontal distances. 
Therefore, the coherence turbulence model shows no signifi- 
cant difference from the one^point spectrum model for this 
landing simulation . 


C. Conclusion 

This study has been concerned with the problems of 
modeling continuous atmospheric turbulence by discrete 
transformation and with application to flight analysis. 
Different spectrum models have been used to generate turbu- 
lence. The results show that the von Karman spectrum more 
accurately describes atmospheric turbulence than does the 
,Dryden spectrum. Both the interlevel coherence of vertical 
layers and the non-Gaussian nature of atmospheric turbu- 
lence, which are important in airplane analysis, can be 
simulated properly by different models. Two discrete 
transformation techniques, FFT and DZT , work fairly well. 
Although the FFT method provides better spectrum simulation, 
the DZT method is faster for computation. 

On the other hand, the disadvantage of the von 
Karman spectrum model is that the non— Gaussian nature of 
atmospheric turbulence cannot be simulated simultaneously. 
For both the von Karman and the Kaimal spectra, only the 
FFT method can be used for the complex irrational filter 


function . 
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Overall, this study has provided a general dis- 
cussion of statistical theories of atmospheric turbulence, 
the detailed description of both random signal analysis 
and digital simulation technique, and the application to 
the simulated turbulence to atmospheric flight analysis. 



BIBLIOGRAPHY 


BIBLIOGRAPHY 


1. Neuman , F., and J. Foster. "Investigation of a 

Digital Automatic Aircraft Landing System in 
Turbulence," National Aeronautics and Space 
Administration TN— D— 6066, Ames Research Center, 
Moffett Field, California, October, 1970. 

2. Reeves, P. M. , R. G. Joppa, and V. M. Ganger. "A 

Non-Gaussian Model of Continuous Atmospheric 
Turbulence for Use in Aircraft Design," National 
Aeronautics and Space Administration CR-2639, 

Ames Research Center, Moffett Field, California, 
January, 1976. 

3. Perlmutter, M. , W. Frost, and G. H. Fichtl. "Three 

Velocity Component, Nonhomogeneous Atmospheric 
Boundary-Layer Turbulence Modeling," AIAA Journal, 
15 (No. 10) :1444-1454, October, 1977. 

4. Frost, W. , and T. H. Moulden . Handbook of Turbulence. 

Vol. 1. New York: Plenum Publishing Corp., 1977. 

5. Hinze, J. 0. Turbulence. New York: McGraw-Hill, 

19 75 . 

6. Kaimal, J. C. "Turbulence Spectra, Length Scales and 

Structure in the Stable Surface Layer," Boundary 
Layer Meterology . Vol. 4. Dordrecht, Holland: 
d 7 Re i del Publishing Co., 1973. Pp. 289-309. 

7. Frost, W. , B. H, Long, and R. E. Turner. "Engineering 

Handbook on the Atmospheric Environmental Guide- 
lines for Use in Wind Turbine Generator 
Development," National Aeronautics and Space 
Administration Technical Paper 1359, Marshall 
Space Flight Center, Huntsville, Alabama, 

December , 1978. 

8. Busch, N. E., and S. E. Larsen. "Spectra of Turbu- 

lence in the Atmospheric Surface Layer, " Danish 
Atomic Energy Commission, Riso Research 
Establishment Report NO. 256, [n.p.] , 1972. 

9. Otnes , R. K. , H. A. Nathans, and L. V. Enochson. 

"A Procedure for Computing Power Spectral Density 
of Gust Data," Report for FDTR Contractor 
MAC 705-03, Wright-Patterson Air Force Base, 

Ohio, September, 1968. 


119 



120 


10. .Mather , G. K. "Some Measurements of Mountain Waves 

and Mountain Wave Turbulence Made Using the 
NAE T-33 Turbulence Research Aircraft," DME/NAE 
Quarterly Bulletin No f 2 , National Research 
Council of Canada, Ottawa, 1967. 

11. Brook, R, R. "A Note on Vertical Coherence of Wind 

Measured in Urban Boundary Layer, " Boundary Layer 
Meteorology . Vol. 9. Dordrecht, Holland: 

D. Reidel Publishing Co., 1975. Pp . 11-20. 

12. Stegen, G. R. , and R. L. Thorpe. "Vertical Coherence 

in the Atmospheric Boundary Layer," National 
Science Foundation GK 35885, [n.p.J, 1978. 

13. Davenport, A. G. "The Spectrum of Horizontal 

Gustiness Near the Ground in the High Winds,” 
Quarterly Journal of Royal Meteorological Society, 
877194-211, 1961. 

14. Chambers, J. M. Computational Method for Data 

Analysis. New York: John Wiley and Sons, Inc., 

15TT. 

15. Hamming, R. W. Digital Filters . Englewood Cliffs, 

New Jersey: Prentice-Hall , Inc., 1977. 

16. Barr, N. M. , Dagfinn Gangaas, and D. R. Schaeffer. 

"Wind Models for Flight Simulation Certification 
of Landing and Approach Guidance and Control 
Systems," Federal Aviation Agency FAA-RD-74-206 , 
System Research and Development Service, 
Washington, D.C. , December, 1974. 

17. Bendat, J. S., and A. G. Piersol. Random Data 

Analysis and Measurement Procedures. New York: 
John Wiley and Sons, Inc . ”, 197l . 

18. Churchill, R. V. Operational Mathematics . New York: 

McGraw-Hill Book Co . , Inc . ~ t 1973 . 

19. Frost, W., and A. M. Shahabi. "A Field Study of Wind 

Over a Simulated Block Building," National 
Aeronautics and Space Administration CR-2804, 
Marshall Space Flight Center, Huntsville, Alabama, 
March, 1977. 

20. Brigham, E. O. The Fast Fourier Transform . 

Englewood Cliffs , New Jersey: Prentice-Hall , Inc., 

1974 . 



121 


21. Pielke, R. A. , and H, A. Panofsky. "Turbulence 

Characteristics Along Several Towers," 

Boundary Layer Meteorology. Vol. 1. Dordrecht, 
Holland: EK Reidel Publishing Co., 1970. 

Pp. 115-130. 

22. Frost, W. , and K. R. Reddy. "Investigation of Air- 

craft Landing in Variable Wind Fields," National 
Aeronautics and Space Aministration CR-3073, 
Marshall Space Flight Center, Huntsville, 

Alabama, December, 1978. 

23. Frost, W, , and W. Crosby, "Investigation of Simulated 

Aircraft Flight through Thunderstorm Outflows," 
National Aeronautics and Space Administration 
CR-3052, Marshall Space Flight Center, Huntsville, 
Alabama, September, 1978. 

24. Frost, W., D. W, Camp, and S. T. Wang. "Wind Shear 

Modeling for Aircraft Hazard Definition," Federal 
Aviation Agency FAA-RD-78-3, System Research and 
Development Service, Washington, D.C., February, 
1978. 

25. Crutcher, H. L. , and L. W. Falls. "Multivariate 

Normality," National Aeronautics and Space 
Administration TN-D-8226, Marshall Space Flight 
Center, Huntsville, Alabama, May, 1976. 

26. Harnett, R. M. "Validation of Automated Uniform 

Random Number Generators for the IBM 7044 and 
the UNIVAC 1108." Report prepared under National 
Aeronautics and Space Administration Contract No. 
NAS8-20082, Marshall Space Flight Center, 
Huntsville, Alabama, July, 1970. 

27. Houbolt, J. C. "Atmospheric Turbulence," AIAA Journal, 

11 (No . 4) :421-437, April, 1973. 

28. Keisler, S. R. , and R. H. Rhyne. "An Assessment of 

Prewhitening in Estimating Power Spectra of 
Atmospheric Turbulence at Long Wavelengths," 
National Aeronautics and Space Administration 
TN-D-828 8 , Langley Research Center, Hampton, 
Virginia, November, 1976. 

29. Etkin, B. Dynamics of Atmospheric Flight . New York: 

John Wiley and Sons, Inc., 1972 . ~ 

. Dryden, H. L. , and A. M. Kuethe. "Turbulence," 
National Advisory Committee for Aeron autics 
Technical Reports. No . 342 . [n.p. , n .n . ) , 1930. 


30 


von Kantian, T. ’’Kantian Theory." Proceedings of 
National Academy of Sciences , U. S . Vol . 347 
{n.p. , n.n. J , 1948 . 

Dryden, H. L. , G. B. Schubauer, W. C. Mick, and 
Skramstad, H. "Turbulence Spectra," National 
Advisory Committee for Aeronautics Technical 
Reports^ No . 581. £n.p. , n.n, ] , 1931 . 



APPENDIXES 



APPENDIX A 


COMPUTER PROGRAM 

A-l, Subroutine GAUSS (IX, SIG, XB, XDP) 

Subroutine GAUSS uses Equation (3-1) , congruential 
method, to generate random numbers. The value of a in 
Equation (3-3) is selected equal to 65539. Also, the 
Central Limit theorem with 12 terms is applied to generate 
Gaussian distributed random numbers. Equation (3-6) . 

Nomenclature 

IX Seed number (input and output) 

-SIG Desired standard deviation (input) 

XB Desired mean value (input) 

XDP Normalized Gaussian random number (output) 

Listing 


SUBROUTINE. GAUSS(IX,5IG,XB,XDP) 
KtAii + 8 XB 
A=0 .0 

DU 30 1 = 1 ,12 • 

IY=IX*6b539 — ) 

IFCIY1 10,20,20 i 

10 IY = 1Y + 2147483647 + 1 RANDU 

20 XP=IY I 

XP=XP*0.46bbol3D-9 i 

1X=I i — 1 

30 A=A+XP 

XDP=(A-6.0)*SIG+XB 

PuTURN 

c.tMD 
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A-2. Subroutine PDF (N, DX, ND, X, P) 

Subroutine PDF computes the probability density 
function of a discrete data set {x^}, i = 1,2,...,N by 
using Equation (4-5) , 

Nomenclature 

N Dimension number for array X (input) 

DX Increment Ax (input) 

ND Dimension number for array P (input) 

X Discrete data set {x^} (input) 

P Calculated PDF of data set {x^} (output) 

Listing 


SUBROUTINE. PDF(li,DX,ND,X,P) 

DlMENSIuN X(N) ,P(ND) ,NP(2048) 
AS=-DX*UND-l)/2 + l ,5) 

DU 30 1=1, ND 
NP(1)=0. 

XS=XS+UX 
XE=XS+DX 
Du 2 0 J = 1 , N 

IF (X C J } . GL. AS. AND . X (d ) . LT . X£) NP(I)=NP(1) +1 
20 CONTINUE 
30 CONTINUE 

DO 40 1=1, ND 
?(!)=!. *NPU)/CDX*N) 

40 CONTINUE 
RETURN 
END 


A-3 . Subroutine FFT (N, NP , NC, XR, XI) 

Subroutine FFT is based on the Fast Fourier Trans- 
form algorithm and computes the Fourier transformation (or 
Inverse) of discrete data. 
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Nome nc 1 a tur e 

N Dimension number (must be a power of 2) (input) 

NP 

NP Power to which 2 is raised, i.e., N = 2 (input) 

NC Control number (input) 

1 for Fourier transform 

2 for Inverse Fourier transform 

XR Real part of the data set (input and output) 

XI Imaginary part of the data set (input and output) 

Listing 


S (J ri P ij U i 1 '4 1- t* K I ( !'i , M P , N C , X R , X I ) 
Dl«r-«SiUw X K ( M , X I C l'i ) 
if l «C. £0.1) GO 10 120 
u u 110 1 = 1 , 0 

no xici)=-xiu) 

12 0 tM l = ;■ j / l 

W 0 1 = I j P - 1 
K = 0 

DU no L=i,NR 
130 DU 1 h0 1=1,02 

R=lHiTK(i\/2**bl)i , UP) 
mKG = o.2 831 LUATU'J) 

C = Cu6 ( ARG) 

U = 5 1 1'; C A R Vj J 
K 1 = K+ 1 
R 1 ii 2 = K 1 + f* 2 

tk=xkcmj^;*c+xicmn2)*s 
T l = AitM.»2)*C-XK(Kl.<2)*S 
XK t K 3 u2 1 =XKtKl ) -1 R 
XilKlri2)=Xi (KI)-TJ. 

XR ( K. 1 )=Xr<Ui ) -UK 
Al CK1 )=XI CM )+'il 
140 K=R+1 
K = R+N 2 

IFCK.UT.R) GO TU 130 

i\ = 0 

WU1=imU1-1 
no W 2 = N 2 / 2 

DU lbO R=l,h 
I=IbirR(K-l ,DP)+1 
if(i.LL.K) GU Tu 160 
1 k=XKCK) 

n=xi(R) 

XR(K)=ARCi) 

XI (K)=X1 C 1 ) 
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XR ( 1 } =TR 
X1(JJ=TI 
160 CONTINUE 

It (NC.LQ.l) GU TO 160 
L>U 170 1 = 1, N 
XRtl)=XKU)/N 
AiWJsXUiJ/rt 
170 CONTINUE 
1B0 KcTUkk 
EiJ O 

FUNCTION 16 XTK(J,NU) 

J1=J 

lulTK=0 

OiJ 200 1 = 1 , NU 

J 2 = J 1 / 2 

l6lXR=ltillR*2 + tJl*-2+J2j 
200 J1=J2 

R L 'l l’> K i J 
El\IL» 


A-4 . Subroutine SPEC (N, NP, DT , XR, XI, SP) 


Subroutine SPEC uses the FFT method to calculate 
estimated spectrum of the input time history, see 
Equation (4-6) . 


Nomenclature 

N Dimension number (input) 

NP Power of 2 (input) 

DT Time increment At (input) 

XR Real part of data array (input) 

XI Imaginary part of data array (input) 

SP Spectrum estimated array (output) 

Listing 


ii li b r(U U 1 1 !'J E ^ t 1 E C C N , l ^ R , D T , X K , X 1 , R ) 
LHi'iEMSlUij aK li. ) , >■ 1 t i* ) , it-’ l 0 i 
h T= C l . * 0'i ) / l-i 
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CALL fr i T C & , rtP , 1 , Af , X 1 ) 
DU 10 1=1, U 

i>F C i ) =KT * ( XR 
10 CUM 11 NHL 
kLTuhh 
EU'u 


A- 5. Subroutine SMOOTH (N, NS, NM, SP) 

Subroutine SMOOTH is the computer program to smooth 
the raw estimate data after discrete transformation. Two 
basic methods, the Frequency Averaging and Segment 
Averaging methods are included. (See Equations (4-9) and 
(4-10) .) 

Nomenclature 

N Dimension number (input) 

NS Number of samples for average (input) 

NM Control number (input) 

1 for the Frequency averaging method 

2 for the Segment averaging method 

SP Raw estimate data (input) 

Smoothed data (output) 


Listing 


sudkoutime smoothcn,ns,nm , SP) 

OlMLi'iSlUN SP(N,NS) 
it'U‘M-2) 100,200,300 
100 w D = N - ii S 

DO 100 1=1, NO 
6UK=0 . 
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DO 140 J=l,NS 
K-l+J” 1 

SuH = 6UM+SP CK , 1 ) 
140 CUN'i' Ii'iUE 

iiPCl , 1 )=SUN/N5 
150 CUilXlwUE. 

GU JO 300 
200 DU 250 1=1, N 
SDK=0 , 

DO 2 40 J = l,t'i3 
6UM=6UM+6PC I ,«J ) 
240 CUNllriUK 

i>PCl , l)=i5UH/WS 
250 CU.'lllrtUb 
300 iibTUKfJ 
biJD 


A-6. Subroutine FILTER (N , NM, DT, Z, V, XR, XI) 

Subroutine FILTER calculates the filter function 
for the different models described in Chapter IV and 
computes Equations (4-13) and (4-14) for the frequency 
components of the output signal. There are three sub- 
routines called by FILTER which are SAL (see A-6 (a) ) , 
COEFl (see A-6 (b) ) , and COEF2 (see A-6(c)). 

Nomenclature 

N Dimension number (input) 

NM Control number for selecting different spectrum 

model (input) 

1 for Dryden spectrum longitudinal component 

2 for Dryden spectrum vertical component 

3 for approximate Kaimal spectrum 

4 for approximate von Karman spectrum 

5 for Kaimal spectrum 

6 for von Karman spectrum longitudinal component 
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7 for von Kantian spectrum vertical component 
DT Time increment At (input) 

Z Height (input) 

V Reference velocity (input) 

XR Real part of the data set (input and output) 

XI Imaginary part of the data set (input and output) 


Listing 


LL dKUUT 1 .'<■ t. F 1 ijTch ( N , NM , DT , Z , V , X R , X I ) 
UirirJMo.hM Xrt(N) ,Xi(L) ,SJ.G(3J , AL C 3 ) 

C ’J !•: .Ml i u / o T l / A 2 , A 2 A , A C 2 , A L» 2 
C i J 1 1 r*.u i / o i / / F w , P 1 , h 2 , C 2 , A 3 
C AL Ij o AL C /j , >0. , vSi G , V ) 

IJ = 1 , / C : m *D'i ) 

L n F = L / 2 + 1 

I F C Nil . Gt . 5 J (jU TO 20 0 

CAUL COct'l (uM,A,C,D,AL»SIG,V) 

DO 100 1=1, uLt 
V = C I - 1 ) + i i* 

F2 = Mf 
A h F = A 2 - F 2 

DD= 1 . / ( A F F -* A F F + A 2 4 * F 2 ) 

H R = D D ♦ C C * A t + A u 2 * F 2 ) 

Hl=uD*(D + Ailr -AC 2 J 

CALL MUTJ CriH,Hl f XRU J#X1(I)) 

100 COimTILUL 
GO l'o 3 0 0 
200 CUwriUUE 

CALL CuLF2(Lrt,A,e,C,AL,SlG,V) 
IFCiiri-6) 270,210,210 
210 Du 2b 0 1=1, NHF 
F=TN*(I-1) 

F 2 = F * t 
KK=A3+f 2 

DD=C (b2 + F2 + C2)**0.S)/(RR**l0.5*Pw) ) 
l'L=ATAiJ.(F/b)+PT*lATAN(F/A) ) 
ilR = DL* (CuS CTrtJ ) 

H I = D D * ( S 1 s 4 ( T H ) ) 

CALL MlJT 1 ( HK , ril , > R ( 1 J , X I ( I ) ) 

2 SO CUriTILUF 
GO TO 300 

270 DU 200 1=1, NHF 
F=('1N*(1-1 ) 

DD=C/(A3+F*F) 

Hrt=uD*A 
H 1 = - D D * F 

CALL MUfl(HR,Hl,XRCI) ,XI(I)) 

2b0 Cunt i lue 
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300 im- = u + 2 

DO 400 1 = 2 , N m F 

J=N1T-I 

XK(J)=XR(I) 

xi cj)=-xicn 

400 COhTINUa 
KC.TUKN 
C.NL) 

ijUbROUTIivt KUT1 (HR , Hi f YR, YI) 
i i rt = »iK*YR-hi *Yi 

U = MK»YH rii*KK 

Y K = Y Y K 
K h 1 U K tl 
t-WD 


A-6(a), Subroutine SAL (Z, AL, SIG, V) 


Subroutine SAL computes the turbulence intensity a 
and the length scale A by using Equations (4-33) and 
(4-34) . 


Nomenclature 

Z Height (input) 

AL Length scale A^, 1 = 1,2,3 (output) 

SIG Turbulence intensity a^, 1 = 1,2,3 (output) 

V Reference velocity (input) 

Listing 


SUbRCJUTI^b SAL(Z,Ab,SIG,V) 

D I M fc. i i & 1 L) N A b C 3 ) , i3 i G C 3 ) 

Z0 = U . 1 

SlG(3)=0.b2*V/ t ALOGtZ/Z0+l ) ) 

SiG(2)=.SlG(3 J/C CO. 5 33 + 0.00 139*ZJ** CO. 8)J 
SlG(l)=SlGt3)/UC.l77 + 0.00 274*Z)**C0.4)) 
IF ( Z.GT . 533 ) GO TU 100 
Ab(n=44„2l + tbuKTC3.2b*Z) ) 

AbCZ J=ALC n 
Ab ( 3 ) -Z 
GU TO 200 
100 Ab ( 1 ) =533 • 

Ab(2)=b33. 

Ab ( 3 ) =533 • 

200 Rc-IuKW 
tt'JL) 
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A-6 (b) . Subroutine C0EF1 (N.M, A, C, D, AL, SIG, V) 

Subroutine COEFl computes the coefficients a, c, 
and d listed in Table 4-2, page 64, for Dryden* s form 
spectrum filter function. 

Nomenclature 

NM Control number (input) 

1 for Dryden spectrum longitudinal component 

2 for Dryden spectrum vertical component 

3 for approximate Kaimal spectrum 

4 for approximate von Karman spectrum 

A Coefficient a (output) 

C Coefficient c (output) 

D Coefficient d (output) 

AL Length scale (input) 

SIG Turbulence intensity o ^ (input) 

V Reference velocity (input) 

Listing 


bUrthuUl lue. CUr.l 1 ( iJ.« , A , C , D , AL , SIG , V ) 
uJLi-lhNSlUU 6 1 G C 3 ) , A L- C 3 ) 

CuKMiji'J /S'l 1 /hi , , AC2 , AD2 

VuA = V/ALC 1 ) 

Pi=3.141b9b 
Pl2 = ra*2. 

Ik- (.oi.hU.A ) GO TU 4 00 
I r t fc - 2 ) 10u f 200,300 
100 A=VOA/P12 

C = SlG(l )/tPl2*PIJ*(VuA**l.b) 

L>=SI6( 1 )/PJL*(VOA**0.b) 

Gu 10 bOO 
200 Vl)A = V /AL C 2 ) 

A=V OA/ PI 2 

C=i*IC;C2)/(FI*PI)*CCVllA/2. 
L>=olG(2)/Pl + ( CVfjA*3./2.J**0.5J 
GO Tu 500 
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300 A=0.155*vuA 

C=0.0432*S1G( 1)*C VUA**1 .5) 
U=61GC1 )*0.279*(VUA**0.5) 
GO TO bOO 
400 M = 0 « 2 fib*vGA 

C=0.1l4*SIGll)*lVOA**1.5) 

D=0.39tf*SiGtl)*(VUA**0.5) 
500 CONTIwUh: 

A 2 = A * A 

A24=4.*A2 

AC2=2,*A*C 

au2=2.*A*D 

tft.TUM 

tNU 


A-6 (c) . Subroutine COEF2 (NM, A, B, C, AL f SIG, V) 

Subroutine COEF2 computes the coefficients a, b, 
and c shown in Equations (4-35) and (4-36) for the von 
Karman and Kaimal spectra. 

Nomenclature 

NM Control number (input) 

5 for Kaimal spectrum 

6 for von Karman spectrum longitudinal component 

7 for von Karman spectrum vertical component 

A Coefficient a (output) 

B Coefficient b (output) 

C Coefficient c (output) 

AL Length scale (input) 

SIG Turbulence intensity ck (input) 

V Reference velocity (input) 
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Listing 


SUoKUUTl.Jt CulFHhK ,A r B,C,AL,SIG,V) 
DiMbi-JSlUM SlG(3),ALU) 

C UN M IJ N / 3 1 2 / P w , PT , b2 , C2 , A3 
V UA=V / AL ( 1 ) 

PI = 3. 141b^t 
Ki2=Pl*2. 

UC:-*M-6) lou,200,300 
100 M = VUA/ (Pi^* 1 . 33S) 

D = 0 . 

C=o.lGCl ) *0.62!>&*( l VUA/P12)**0.333 333) 
P n = b . / 6 . 

G U TO 400 

200 AA=0. Ib4*( Cl./0.041)**(b./3.)) 

A6=S0kT l A A ) 
bbsSliRKO. 164/0.041 ) 
m=(VUA* + (5./6. ) )/AS 
b = 0 ♦ 

C=SiG( 1 )/A3*(V0A**0.33333)*BB 
Pw=5./6. 

GU TU 400 
300 VuA=V/AL(2) 

Pw = l 1 . /b. 

A=VUA/(P12*1 .339) 

B=3URT(2.*SIG(2)**2/V0A*(A**tPW*2) ) ) 
C=B/(0.6123*A) 

400 PT=-PW 
A3=A*A 
62=b*B 
C2=C*C 
Rtl'UKN 
tM> 


A-7. Subroutine DZT (N, NM, DT, Z, V, X, Y) 

Subroutine DZT computes Equation (4-15) for 
different models, and this is the Z— transf ormation tech- 
nique for the filter system. 
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Nomenc 1 a tur e 

N Dimension number (input) 

NM Control number for setting different spectrum model 

1 for Dryden spectrum longitudinal component 

2 for Dryden spectrum vertical component 

3 for approximate Kaimal spectrum 

4 for approximate von Karman spectrum 

DT Time increment (input) 

2 Height (input) 

V Reference velocity (input) 

X Input array 

Y Output array 

Listing 


SUBR0UT1 im£ DZT(N,HM,DT,Z, V , X ,Y) 

LlriblvJSIOu X(NJ , Y(N) ,S1G(3) , AL(3) 

PI = 3 , 1 4 1 b93 

CALL SALCZ, AL ,LIC, V ) 

CALL COct'l (Lrt,A,C,L, AL,S1G, V) 
Cl-2.*bArH-2.*PI*A*DT) 

C2 = -(fc.XP(-4.*PI*A*D , n ) 

b=LO«T ( 2 . *P1 ) 

Cl 2*Cl/2. 

CbA = C/ (ci*A*A ) 

ACA=( /A 

L)l=L>+(CbA + Cl2*(ACA*0T“CBA) J 
D2 = L'i*Cl^^CCbA»tC12-l • )"ACA*L)T) 

Y( 1 )=0. 

Y (2)=0. 

Uu 30 1 = 3, N 

Y(JU=Cl*YCl-n+C2*YU-2)+f)l*XU-l) + D2*A(l-2) 
30 CUiiTXiiUt 
RETURN 
E.N0 
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A- 8 . 


Subroutine INPUT (N, NM, DT, IX, XR, XI) 


Subroutine INPUT generates the random number array 
{x^}, i = 1,2, . ,.,N, which standard deviation a is adjusted 
for FFT and DZT, as shown in Equations (4-25) and (4-32) . 


Nomenc lature 


N 

NM 


DT 

IX 

XR 

XI 


Dimension number (input) 

Control number for selecting different method 
(input) 

1 for FFT method 

2 for DZT method 

Time increment At (input) 

Seed number for calling GAUSS (input and output) 
Real part of the array (output) 

Imaginary part of the array {x^} (output) 


Listing 


SUoKUUTlNt INPUT(N,NM,DT,IX,XR,XI) 
D X ft c i'i 1 u “J X P C N ) , X I C N ) 

IKNrt-2) 100,200,300 
100 iiIL = SQRT(0.b/DT) 

00 TO 300 

200 SiO=SwkTC3.l41b56/L)T) 

300 DO 400 1 = 1, N 

CALL GAUSS(IX r SlG,0. ,XDP) 

Xrt C 1 ) = XLP 
Xl(l )=U. 

400 CukllNut 
Kt. l Oku 
LNO 



137 


A- 9 . Subroutine NONGAU (N, NC, DT, Z, V, R, XR, XI) 

Subroutine NONGAU generates the non-Gaussian turbu- 
lence based on the model of Reeves, et al. [2], and 
calculates the filter functions listed in Table 2-1, page 
26. There are three subroutines called by NONGAU, which 
are INPUT (see A8) , C0EF3 (see A-9(a)), and FFT (see A-3) . 

Nomenc lature 

N Dimension number (input) 

NC Control number for specifying the velocity 

component (input) 

1 for longitudinal component 

2 for lateral component 

3 for vertical component 

DT Time increment (input) 

Z Height (input) 

V Reference velocity 

R PDF adjusting coefficient r (input) 

XT Real part of simulated turbulence (output) 

XI Imaginary part of simulated turbulence (output) 
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Listing 


ii U rt K u U 'll h L NQNGAU(N,NC,DT,Z, V,XR,Xi,R) 
DIMENSION XR(N) , X I C K ) ,AlA( 3J ,51GS(3) ,UK(2Q48) 
CuMtfUi-l/6I7AL ,51G, IP, A2 , A24 , ADZ , AC2 
tit 2= ( AhUG ( 1 . +h ) ) / ( ALOG C 2 . 3 ) 
roi = fi/2 

Ti'i = 1 . / ( N *UT ) 

WhF=ijH + 1 

call salcz, ala,sigs, V) 

AL=ALA(WC1 
61G = £>1G6 ( iN C ) 

UK=1 ./(S-JHTt 1 ,+R*hJ ) 

KUK = R*Lik 
NK = 1 

100 CALL COEr 3(Z, V,A,C, D,NC,MR) 

CALL 1 N PUT ( N t 1 till' , 1 X , X K , X I ) 

CALL FFTIN,N12,1,XR,XI) 

DU 150 X = 1 # tx HF 
F=CI-1 )*TN 
F2=F*F 
AMF=A2-F2 

UU=1./(AMF*AMF+A24*F2) 

hK=Up*(C*AME+AD2*F2) 

Hi=DD* (D*AMF-AC2 ) *F 
XXR=XR(l)*hfi-Xltl)*Hi 
XI(l)=XKll)*hI+Xl(I)*HR 
150 XRCI)=XXR 
WTT=N+2 
DO 190 1=2 , Nh 
J=NTT-I 
XK(J1— XRC1) 

190 X1(J)=-X1(1) 

CALL FFT Ci'< , NT2 , 2 , XR , XI ) 

IFCaiK-2) 200,300,400 
200 DU 250 1=1, N 
250 QRU) = XR(I) 

NR = jNR + 1 
GU TO 100 
300 DO 350 1=1, N 
350 DR ( 1 )=URll)+XR(l) 

Wk = 3 

GO TU 100 
400 DU 450 1 = 1, N 
450 Xr(1)=XR(I)*UK+QR(I)*KOR 
Ke-TU RN 
ELD 


A- 9 (a) , Subroutine COEF3 (Z, V, A, C, D) 

Subroutine COEF3 computes the coefficients for the 
filter functions listed in Table 2-1, page 26. 



Nomen c 1 ature 


2 Height (input) 

V Reference velocity (input) 

A Coefficient a (output) 

C Coefficient c (output) 

D Coefficient d (output) 

Listing 


6UdK0UTli(fc: cutt'3 tz, V , A,C,D,NC,iiH) 

G U rt rt iJ u / GT / A L. , £> i G / 1 P , A 2 , A 2 4 » A D 2 , A C 2 
VUA= WrtL> 
w=b:)KT(. 1 2b) 

VUA2=0.5*VuA 
U'(uc-2} 100,200,200 

100 1KUK-2) 110,140,170 

110 A= \ZGA2 

C=4.*SIG/V0A*A*A 

0 = 1 . 

GU TU 400 
140 A=VUA2 
C = A*A 
0=1 . 

GO TO 400 
170 A= ULi 

C=SIG*tSGKT(2./VuA))*A*A 

0 = 1 . 

GO To 400 

2 00 lt(.Uii-2) 210,240,270 
210 A = V 0 A 2 

C=olG*rt/lVuA**2)*A*A 
0=1 . 

GO To 4o0 
240 A=VUA2 

C-U. 

0 = A 

GO TU 400 
270 A= Vu A 

C=S1G*A/ (Soul C Vt*A ) ) 

D = G* t . 732/ VuA 
400 A 2 = A * A 

A2 4 = 4 . * A 2 
AC2=2 . + A*C 
A I) 2 = 2 . * A *0 
kr/i'URN 
c,Ni) 
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A- 10 . Subroutine COHER (N, Z, DZ, DT, Yl, Y2) 

Subroutine COHER simulates two levels of turbulence 
Yl and Y2 with the coherence relation described in Equation 
(4-43) . There are five subroutines called by COHER: 

BETA (see A-10 (a) ) , NTRA (see A-10 (b) ) , VHAT (see A-lO(c), 
FILTER (see A-6) , and FFT (see A-3) . 

Nomenclature 

N Dimension number (input) 

Z Height of the level one (input) 

DZ Vertical separation for two levels (input) 

DT Time increment (input) 

Yl Simulated coherent turbulence at level 1 (output) 

Y2 Simulated coherent turbulence at level 2 (output) 


Listing 


SUokUUT i CUHEk(N f Z,LZ,L)T,Yl » Y 2 ) 

01rtr,koJL0rt Y1UO , Y2UO 

DlrfLNSluu K ; i K ( 1 1 , b 1 2 ) , R ivl 1 1 1 1 » b 1 2 ) ,Xk(bl2J » X I ( b j 2 ) ^k(ll ) 

Cui'iMu,\l/ r2/t>i ,Z1 , DX , LPO , A 

C u .'i in U R / X 3 / l'i P 2 1 , k H F , L T 2 

UPt,r=5 . 0 

IX=6bb49 

ui 2 = ( ALUG C 1 . +;» ) ) / C ALUG C 2.) ) 

ntl — ii/ 2 

i'i H e L H + 1 

im P 2 1 = 1 1 

Pl = 3. 141593 

A=17./(2.*PI) 

Ll-L 
D Z H = 3 , 

Ib=l 

tPU=l')ZM*Pl/(UKbK*L>T) 

L)A = UkbF*0T/ C P 1*2 . ) 

C Abb dETACf.Pil , B) 

CALL NTKA(N,l)T,IX,RHRfRNI,XR,Xl) 

100 CALL VHAT(N,DT,RilK,RNl,XR,Xl,n) 
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CAut, F I h T t K ( N , U c \ , 0 T , Z , U K h F , X K , X 1 ) 
CALL F FT(il,UT2,2,XR,XIl 
IFC1L.E0.2J GU TO 300 
UU 200 K=1 # H 
tflCKJsXKCKj 
200 CUNTiuUfi 
4=4-i)Z 
iL = iL+l 
GU i'U 100 
300 UU 310 K = 1,N 
i 2(t)=XK(K) 

310 Cuul'lwUE 
8LTUK.M 
fc.ND 


A-10(a). Subroutine BETA (N, B) 


Subroutine BETA computes the coefficients B^ of 
Equation (4-38) using Equations (2-45) and (2-42) . 


Nomenclature 

N Dimension number (input) 

B Coefficient array B m (output) 

Listing 


SUoKtjUTINe, 8LTAUJF21 ,b) 

0 i p | 5 1, i * o 1 U i'J b (.uP 21) 

CurfMWT2/Pl ,21 ,l)X,LPU, A 
Ut»l = OP21-l J/2 + 1 
U=A*tPu/2. 
dCNPl 1=0. b 

6Ul-i= 1 . 

LXhUstXP ( -L>) 

Pbz=8l/D 

ANU5isCl.-t-c.AL0J/ U ,-LXLOJ 
UL 300 l=l,uPt,2 
K=NP1-1 

uCKJ=ANiJH/Cl .+(8t2*ll**2l 

IFCK.LQ.U bu TCI 100 

b(K-l )= l./Cl.+(Pt2*CI+l))**2) 

SUM=SUKi + 2.^(tiCK) + B(K-l ) 1 

b(K-] Js.SOKTCeiCK-l 1 1 

GO 10 200 

100 DU i v i = b U A + 2 . *d ( K) 

200 o(M=vsUKTCbUM ) 

300 CUiMliLUt 
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l»u 400 I - 1 / N P 1 
400 b ( 1 ) = b ( I ) X £>uH X ( iSli tf ) 
ht.TURI') 
lNO 


A-10(b). Subroutine NTRA (N, DT, IX, RNR, RNI, XR, IX) 

Subroutine NTRA generates 11 sets of random signals 
for the coherence model. 

Nomenclature 

M Dim nsion number (input) 

DT Time increment (input) 

IX Seed number for calling GAUSS (input and output) 

RNR Real part of the random signal (output) 

RNI Imaginary part of the random signals (output) 

XR, XI Working space 


Listing 


SUdRUUTIhl. iVfKA(h,DT,IX,RNR f RNl,XR,Xl) 

OlKiL.Jb 1U.N RljRtll / b 1 2 ) , hfl i ( 1 1 » b 1 2 ) , XK 1 b 1 2 ) / XI (512) 

CUi-iKUU / T 3 / N P 2 1 , M H f , N T 2 

DU 300 J=1,LP21 

CALL INPUT Ui, 1 ,DT,iX,XR,XI) 

CALL ttT(L,N‘l2, 1 ,XR,XI) 

DU 200 1=1, NHF 
RMRUl, I)=XR(1 ) 

RHl(J,I)=Xl(l) 

200 CUuT 1 uUfc 
300 CuNI-InLL 
RETURN 
ti\i D 
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A-10(c). Subroutine VHAT (N, DT, RNR, RNI , XR, XI, B) 

Subroutine VHAT is a part of the coherence model, 
and computes the X* (k,z) in Equation (4-39) . 

Nomenclature 

N Dimension number (input) 

DT Time increment (input) 

RNR, Real and imaginary part of random signal (input) 

RNI 

XR, XI Real and imaginary part of calculated X' (output) 

B Coefficient array in Equation (4-38) (input) 

Listing 


oUnHuUTXNE VnAT (M , UT , RNR , RW I ,Xfc,Xl,b) 

lUrtLwSlUrt RhKi 11 ,512) ,KrU(U ,S12) ,XkCbl2) ,X1(512) ,u(ll) 

CUrtri0jn/T2/Hl , L 1 ,UX,i.Pu # A 

COnttUrt/l'i/hf’il ,vU:F ,MT2 

rtl 2=uP21+l 

Pl2=Fl*2. 

CK1T = J1X*h*P1/10. 

A k. 1’ A = P I * 2 . / C JLPG*DX*N ) 

DU 200 K = 1 , N H F 
FKK=1 ,*C«W ) 

lF(t KK.Gb.CklT) t-kk = CKlT 

xhckj=o. 

XI CK)=0. 

DO 100 h— 1 , i\ Pl 
AUG = At.T A* Q wP 1 - Ij ) 

AUGfc=(-AuG) KK 
rtUGF—CAliGJ-rt'KK 
C6N = C0S C AUGhf ) 

C6P = C0* C AuGF ) 
bi»N = SlM C AuGU ) 

5NP = Sltt t AJuP) 

XK(K)=XRlft)'t(Kl'JK(Nl2~aj,K)*C. , £>ft'tkMt<(Lf,K) +C6P-K N 1 C H 1 2-L , K ) * 
♦ iiWrJ-K,MlCl>,K) *4>wP) *b(L) 

Xi(k J=Xi (n) + (kkl ([112-L,K)*Ci>N + khil t L , K) ♦CbF + hf'iK ( M 2-b , K ) ♦ 
+ 6k t* 4 Kim K C G , X ) * 61 * F ) * b ( L> ) 

100 CUh'IINUE 

200 C'J.Vl 1 riU t, 

Ktl’t'KI-J 
t,H L) 
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A— 11 . Subroutine COHERS (N, NS, DT, Y1 , Y2 , CO) 

Subroutine COHERS computes the coherence function 
of the two time histories Y1 and Y2 by using Equation 
(4-41) . 

Nomenclature 


N 

Dimension number (input) 


NS 

Sample size (input) 


DT 

Time increment (input) 


Yl 

Time history at level 1 

( input) 

Y2 

Time history at level 2 

(input) 

CO 

Estimated coherence of Yl and Y2 (output) 


Listing 


oUdkUUTI KC, CUHbKb ( w , NS , DT , 1 1 , i i , CO J 
Ol'-'LluSlUi'. Yl(n.S r L) , Y2 (.;mS,i'] ) , CU C IJ J 

UlKhi.oIuU Xk C 204 8 ) , XI L 20 A 8 ) , CuK ( 102*) , CUI t 1 0 2 4 ) , SP tl 0 2 4 ) 
D i ft l ri 6 1 0 ii CukK(10?4) ,CUIi(1024) 
ivX2 - ( ALU<j C 1 . *1\) ) /AijIJGC 2 . ) ) 

Nri = N / 2 
Kr = uX/l'iH 
UU y0 K=1,.NH 
CORK C K ) =0 , 

Cull (i\)=0. 

90 CUriXIHUt 

OU 3 00 L=1,NS 
i)u 100 K=1 ,ti 
AftCh )=Y1 CK f L) 

XI (K )=0. 

100 CuuTlNUE 

CALL SPEC (N , NT2 , LT , XR , X I , SP ) 

LU 110 K=l,Nri 
CUK t K ) =XK CM 
COI IK) = Xl t K ) 

¥ 1 ( K , L ) =SP ( K ) 

no co.m.vUK 

LU 200 K~ 1 , N 
XktM = Y2tA,L) 

XI ( K ) =0 . 

200 CUnTIiMUE 

CALL SPEC(N,NT2,DT,XK,XI , SP) 

L>U 210 K = 1 , W H 
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CuKR(K)=CUKk(K)+CuR [K)*XK(K;-CQ1(K)*XIIK) 
CulI(K}=CuiI(K)+COR(KJ*XI(K)+CLil{.ls)*XKtKJ 
Y2(K,L)=SR(K) 

210 CUHTIKUE 
300 continue 

L>U 4 0 0 K=1,NH 
SUK1 =0 . 

SUK2=0. 

00 3 bO b=l,NS 
SU,ii=SUMl + i 1 (.K,L) 

£>UM2=SUM2 + Y2tK,L) 

350 CUNTXiaJK 

1 1 (K, 1 J=oUhl /»\S 

1 2 CK , 1 ) =SUH2/Uii 
CiJKh(K)=CUHk(iO /foS 
CUIl(KJ=CUj.i(N)/Nb 

40 0 CUM IN UK 

OU 500 K=l,kH 

CuCM=KT*rtr4 ICURF(K) **2 + C011UO**2)/(U (K» l)**2(h, 11 ) 
50 0 CUNilNliL 
k Li' U Km 
C- iM L> 


A- 12 . Subroutine WIND (XP, ZP, T, V, KCK) 


Subroutine WIND computes the turbulent wind velocity 
and U 2 (see Equations (5-1) and (5-2)). 


Nomenclature 

XP , ZP Aircraft position (input) 

T Time (input) 

V Aircraft speed (input) 

KCK Starting number KCK = 1 (input) 

COMMON/WINDS/ 

WX ,WZ Wind velocities U^, U 2 (output) 

WXT, WZT Time gradient terms SU^/St, 3U 2 /3t (output) 

WXX , WXZ Spatial gradient terms 3U-^/3x, 3U^/3z (output) 
WZX , WZZ Spatial gradient terms 3U 2 /3x, 3U 2 /3z (output) 
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Listing 


SUBtUJlITlNt; WIND (X?,ZP,T,tf ,KCK) 

CUriMiji./ia:iL-3/WX,WZ,tiXT,WZT,vfAA,WXZ,UZA,WZ£,ZC>,U.>rAA 
COMmiN/STSK/FAT.LDM I 

31 FOR >i Al" (2 a, 1 ***PUJ.iT X = ' ,t,13.h,* Z = * ,£13.6,' 1$ OUTSIDE ' ) 

AKP=0.4 

1FCXP.LI.0. .OK.ZP.LT.O, ) GQ TO 21*0 
UUK = ’./DRl*Ui»‘I Ah/AK^' ' 
l»X--uOKf (AlOG( (ZP + Z.J) /Zli) ) 
i« Z = 0 . 

*« A A — li . 

wXZ = -i' IK/lZP+Zb) 

WZX=0. 
wZZ=0 . 
tfXT = 0 . 

WZT=0. 

CALL. VLLO(XP,ZP,T, V ,KCK) 

CO TO 215 

210 wRiTLC 6/31) XP,ZP 
215 HfcbTUrtN 
END 

SUBROUTINE VbLOCXP, ZP,1, V.KCK) 

Cum ! ,i i >1 / T1 1 / OX , OT , k« r. , k T 2 
COMMON / XT 2 / GG ( 2 , 1 2 8 ) 

COMMON /XT 4 / 1 X , 8 R , d ET H 
COMMON/ TT 5 / XPK , ZPR , T END, ZOWR 

CQrtm)ll/XI..DS/rfX , MZ , riXT , WZT, Vi XX , V XZ , w z X , t.ZL, Zi.>,V6 fAK 
ONii— 

WXk=A8S(«X) 

IFCABSCWX) .LT. 0.0000001) WX 8 =0.00001 
IF(KCK-l) 101,101,103 
101 XPR=XP 
ZPR=ZP 

CALL RUMDUM C ZP ) 

Zm/R = ZP/WX8 

Tbf.u = iXiN 

GO iu 105 

103 ZOw=ZP/WXR 

If (ZP.LT.20.) GO TO 105 

RRR=ZUW/ZDivK 

DDD=AriS ( 1 -RRR ) 

1FCDD0.LT. 0.2) GO TO 105 

104 CALL KUMDUHCZP) 

ZOWR=ZOW 

XPR=XP 

ZPrt=ZP 
TENO=T + l)NN 

105 TKP=T+ABS( (XP-XPR)/V) 

1F(TPP-Tc..jD) 109,104,104 

109 tR=C THP+U, mN-TEHD)/DT 
NT=ERf 1 
NTP.=rtT* 1 

If CNTP.GT. 128) GlJ TU 104 
UP=GG( 1 ,nT) 

VP = GG(2,.ir) 

A'X=WX + uP 
wZ=WZ* ^P 
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dld=v*d r 

WXX=V/XA ♦ CGGl 1 ,NTr*)-GG( l.ND )/0L0 
\-)‘/..K=wZXi (G6(?,NTP) -GG( >,Mn ) t OLD 
1 2o C0>-i1 I NO £ 

RETURN 

fc.NO 

auiiKuijTiNfc: hum dun i 7,p ) 

CUiiMi’iN/TT 1 /UX , DT , NNO , «T2 
CU.iNtt,J/TT.:/'C.U( 2, l r- j 

Cu.m i-tijij / T 1 4 / 1 x , R it , t t v h ■ 

CJK.-UWFT/XGKl 12b ) , X G I ( 128) 

Curt AOfl/STSU/t'AT , Vi D H .1 

COMMON/ WluDS/WX , WZ , *iXT , iflZT , V. XX , X Z , »ZX f WZZ , Zi i , U3TAF ' 
SIvjX =;>nKT ( J . I4l59b/DT) 

OU 110 1=1,2 
DO 105 0=1, M N 1 i 
CAUL GAU35t XX , SIGX ,0. ,K) 
aGK 1 0 ) =k 
J.0 5 CU.jT1.JUc. 

UU 103 J = 1 , N N N 
XGI (J) =0 . 

108 CONTINUE 

CALL F f T ( N N N , N T 2 , 1 ) ' 

CALL F XLTcK ( ZP , I ) 

CALL FFT(NHN,KT2. 2) 
wEITSCo,!) CXGRCK) , K= 1 , NNM , 5 ) 

1 FORMAT (IX, 20F5 . 2 ) 

DU 109 J=1,NNN 
GGC1 , J ) =XGP. ( J ) *FAT 

109 COfjT IiiUE 

110 CON Tin US 
RETURN 
END 

oUBrtOUTluE FILTERCZP, IOP) 

COMMON / TT 1 / DX , DT . N »'J ri , NT 2 
CUMM0N/FT/XGKC128) ,XGI (128) 

COMMON /W IN DS/WX , WZ,WXT,WZT,WXX,VjXZ,WZX,wZZ,ZO,USTAR 
PI=3. 14159b 

AL = ZP/(0. 177+0. 832/304. 79*ZP)**1 .2 
UMEAN=ABS(WX) 

IFCUMEAN.LT. 0.0001) UMEANsO . 000 1 
5IG3=0 . 52 ♦ui’iEAL / i ( ALuG(ZP/ZO + l ) ) ) 

SlG=SIG3/t (0.177+0.00274*ZP)**(0.4) ) 

IFClUP.EO.2) SIG=S1G3/C (0.58 3+0. 00 139 *ZP)**( 0.8) ) 

Al)V = AL/UMEAN 

MH-NiJN/2 

NHF=NH+1 

TN=1 ./CNNN*DT) 

100 CONTINUE 

1FC10P.E0.2J GO Tfl 115 
Al=l./(1.339* aOV * 2 . * t> 1 ) 

Cl=4.*At)V*blG*£»lG*lAl* + l.oooooo) 

PO=5 . / 1 2 . 

C ! S-50RT (Cl) 

F AC=-5./6. 

f\ l 2 = A l * A l 
DO 110 K= 1 , HHp 
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S 2 = o S 

l) = ClS/t t A t 2 -» £> 2 J + + t-'O ) 

A N G = fr‘ A C * ( A T A f* ( r> / A L ) ) 

HksD* (COSCAJ'G) ) 
hi=D*(Si;j(At;Cj ) 

XAK-XGk(XJ ♦HR-XG I (K) *H I 
XG1(KJ=XG1(K) M-m< + av.;kIK ) * H 1 
AGkU()=XAR 

no continue 

go 10 300 

i 1 5 CONTINUE 

A=1./(1.339*A0V*2.*R1) 

C = A*(*S(jRTC3./6. ) ) 

C2=C*C 

b2 = 2.*SlG*SlG*AUV*(A**Ul ./3 . ) )/ iCi ) 
6=3 QRTC62) 

A2=A*A 
f h = i 1 . / 1 2 . 

Ra2 = 1 1 ,/6. 

00 120 K=l,riHF 
S= C K- 1 ) *TN 
S 2 ^ S ^ S 

D=d*CSOKT(C2 + S2) )/( ( A2+S2)**iM) 

A hi G = A T A N C 5 / C ) - P W 2 * C A T A N C S / A ) ) 
rtR=D*(CUo(ANG) ) 

Hl=D*(SI»(ArjG) ) 

XAA=XGR(n)*hR-XGI (K)* ttl 
XGl(rO=XGI(.Kj*hk + XGH(K) ’■hi 
XGR C K) =XAA 
120 CONTINUE 
300 H3=Nfirf+2 

DO 330 J~2 » NH 
K=w3-J 

XGRCK)=XGRCJ) 

XG1 ( K) =-XGI C J) 

330 CONTINUE 
KETURu 
fchD 



APPENDIX B 


VARIANCE OF DZT 

An approach to derive exact equations for the out- 
put variance by using the DZT method is to calculate the 
variance directly by expressing the difference equation 
involving input and output variables in terms of the input 
variables only. Recall Equation (3-46) : 


y n+l = Ay n + Bx n 


where 


A = exp(-aAt) 


B = 


/B 


a = -r 


c = 


2 Vo 

~ 7wT 


a 

2 


( 1 - exp (-aAt) ) 


Rewrite Equation (B-l) so that 


n+1 


= Ay + Bx 


n 


n 


y = Ay . + Bx , 
2 n Jr n-1 n-1 


y 2 = Ay x + Bx L 

y 1 = 0 . 


(B-l) 


(B-2) 
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Substitute into y Y 2 i n to y^,-.., etc., which results 
in 


Y n+ 1 = A n 1 + A n " 2 Bx 2 + 


+ Bx 


n 


n 


= I 


i=l 


a . x . 

1 1 


r 


where 


a . - B 
1 


i-1 


Since x^ and y^ are independent, the variance can be calcu- 
lated by 


n 


a = 


a lim J af 
x n+°° i=l 1 


.2 c 


-At 


n 


07 — -jr ( 1 — e ) lim £ ( e~ aAt ) 


i-1 


x 2 
a 


n^°° i=l 


(B-4 ) 


Note that 


OO 

. L b - n 


i=0 


Therefore, Equation (B-4) becomes 


o ~ 2 

a 2 

X 


2a y (l - e - aAt ) 2 

tt a ” "-2a At 

1 - e 


(B-5) 



Consider the case for aAt very small 
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2 ^ 

( i _ e _aAt ) _ ( 1 - 1 + aAt - + . . . ) 

! _ e ~2aAt 1 - X + 2 aAt - . . . 

(B-6) 

aAt 

2 * 


Hence , Equation (B— 5) reduces to 


a 2 

a 

_Z 

2 

a 


x 


At 

TT 


a 


2 

y 


(B-7) 



APPENDIX C 


EQUATIONS OF MOTION FOR LANDING SIMULATION 

The three-degree-of -freedom model for aircraft 
motion presented in this appendix follows the general form 
developed by Neuman and Foster [1] , except the non-linear 
terms are retained. It accounts for both vertical and 
horizontal mean wind components having both time and 
spatial variations . 

Figure C-l illustrates the forces acting on the 
aircraft. These include: 

F^ = thrust of the engines 

L = lift 

->■ _ 

D = drag 

U = wind vector 
mg = gravitational force . 

Figure C-l shows the orientation of the forces with 
respect to the velocity relative to the earth V, the 
velocity relative to the air mass V , and the fuselage 
reference line FRT of the aircraft. 


152 



*! 
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Figure C-l. The forces acting on the aircraft. 

The basic variables are V, y, q, a', x, and z, 

where 

y = angle between V and x-axis (flight path angle) 

q = time derivative of 8 (pitch rate) 

a'= angle between V and the FRL (angle of attack) 

a 

From a direct force balance along each direction 
and basic kinematic relations, the equations of motion can 


be derived as 
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V = -Dj ( oos 6 + c sin 6 ) sin y + S ^ J1 Y + D g F T 003 ^ + a ^ 


U-i V 

1 a 
V 

(C L 

— D F m 
7 T 

+ D 5 


( c cos 6 - c D sin 6 ) - D 2 


. 99 s T + D F 

V 6T 


sin ( 6 T + a ) 
V 


- 1- W ,« 2 ^ - Vi - f (ft .in V .4, =»,' > 


x = V oos y 
z = - V sin y , 


where D^, D 2 , . . . are the nondimensional coefficients 
and c^, c^, and c^ are the aerodynamic coefficients defined 
by 


D, 


= C L +C L a ' +C L, 6 E + V- lG L q+C L a ' ]+C L 


o a 


a q 


GE 


C D = C D +C D +C D2 a ' 2 + C D rr 
o a a GE 


D, 


c M= c M + % S E + V~ l \ * + \ a ' ] + 

o a 0^ a G[ a Til 

hi 


The wind is seen to enter the equations in the form 
of a gradient or wind shear and U 2 . The expanded form 
of these equations is: 
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3U 1 + 

9U 1 

dx , 

9U 1 

dz 

TF + 

9x 

at + 

9z 

dt 

au 2 + 

3U 2 

dx 

9U 2 

dz 

9t + 

9x 

dt 

9 2 

dt * 


Thus, both spatial variations and temporal variations in 
atmospheric motion influence the equations in the wind 
coordinate system. 


APPENDIX D 


LANDING RESULTS 

This appendix presents landing results for non- 
Gaussian, Dryden spectrum turbulence, and Gaussian, von 
Karman spectrum turbulence (Table D-l) . 
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Table D-l. Landing Results 


No. 

Non-Gauss ian , Dryden 
Spectrum Turbulence 

Gaussian, von Harman 
Spectrum Turbulence 

Landing 

Position 

Standard 

Deviation 

Landing 

Position 

Standard 

Deviation 

1 

7033.13 

62.41 

6926.15 

-35.97 

2 

6924.73 

-46.0 

6984.11 

21.99 

a 

6983.23 

12.51 

6991.46 

29.34 

4 

6933.35 

-37.37 

7009.88 

47.76 

5 

6968.89 

-1.83 

6986.71 

24.59 

6 

6998.80 

28.08 

6992.73 

30.61 

7 

6977.84 

7.12 

7022.72 

60.60 

8 

6917.63 

-53.09 

7016.62 

54.50 

9 

7007.27 

36.55 

6909.70 

-52.42 

10 

6920.63 

-50.09 

6923.54 

-38.58 

11 

6976.52 

5.8 

6962.78 

0.66 

12 

6939.20 

-31.50 

7008.04 

45.92 

13 

7057.49 

86.77 

6951.11 

-11.01 

14 

6974 

5.03 

6948.43 

6.31 

15 

6969.15 

-1.57 

6914.48 

-47.64 

16 

7009.27 

38 . 55 

6925.56 

-36.56 

17 

6994.06 

23.24 

6927.14 

-34.98 

18 

6944.93 

-25.79 

6896.60 

-65.52 

19 

6946.49 

-24.23 

696^. 51 

-1.61 

20 

6936.05 

-34.-67 

6964.12 

2.00 

Average 

6970.72 

37.51 

6962.12 

32.43 


: u + = 0.5. z = 

71 O 


Note 


0..1 
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